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1  Introduction 


Filtering,  prediction,  and  smoothing  (FPS)  are  the  three  basic  components  of  the  data 
assimilation  process  in  target  tracking.  An  analytical  solution  of  the  FPS  problem  is  possible 
only  in  a  handful  of  particular  cases,  the  most  important  of  which  is  linear.  In  this  case 
the  solution  is  given  by  the  Kalman  filter.  However,  in  many  important  cases,  such  as 
passive  sonar,  radar  warning  systems,  infrared  search  and  track,  the  systems  are  generically 
nonlinear.  To  date,  the  extended  Kalman  filter  (EKF)  has  been  the  dominant  algorithm 
technology  in  real-time  estimation,  tracking,  and  similar  applications.  A  major  reason  for 
its  success  has  been  the  fact  that  it  has  offered  a  reasonable  compromise  between  real-time 
operation  and  satisfactory  performance  in  some  nonlinear  problems.  On  the  other  hand,  the 
EKF  is  a  completely  heuristic  algorithm,  requires  readjustment  to  each  particular  problem, 
and  is  unstable  in  nonlinear  problems  which  involve  jumps,  maneuvers,  etc. 

Nonlinear  filtering  is  the  process  of  computing  estimates  of  the  current  state  Xt  of  a 
nonlinear  dynamic  system  (e.g.,  a  rapidly  maneuvering  target),  given  current  measurements 
(which  are  nonlinear  functions  of  the  system  state)  together  with  some  degree  of  knowledge 
of  the  states  of  the  system  at  previous  instants.  The  state  may  include  such  unknown  target 
characteristics  as  position,  speed,  acceleration,  aspect  angle,  etc.  Kinematic  data  Zk  are 
collected  at  discrete  time  instants  k  —  0, 1,  2, . . ..  The  relationship  between  measurements 
and  target  states  is  modeled  by  a  (generally  nonlinear)  measurement  equation  of  the  form 
Zk  =  h{xk,  Vk)  where  Vk  denotes  the  noise  process.  The  expected  range  of  possible  behaviors 
of  the  target  is  modeled  by  Markov-state  transition  equations  of  the  form  =  b(xt,Wt) 
where  Wt  is  another  noise  process. 

In  many  practical  situations  the  EKF  and  its  variants  do  not  perform  well  when  the  model 
functions  b  and  h  are  highly  nonlinear.  The  reason  is  that  the  EKF-like  methods  approximate 
the  true  posterior  probability  densities  Pk\k{^k  |  Z^)  of  the  target  state  by  Gaussian  laws 
(Z^  =  {zi^ . . .  ,Zk)  being  the  data  observed  up  to  time  k).  When  a  target  executes  a 
maneuver,  Pk\k{^k  |  Z^)  often  becomes  multi-modal,  with  different  modes  corresponding  to 
the  most  likely  time-varying  alternatives  of  the  current  target  state.  As  data  are  collected 
one  of  these  modes  will  eventually  dominate  the  others,  corresponding  to  the  actual  state  of 
the  target.  By  contrast,  EKF  is  forced  to  approximate  pfc|jt(a;jfe  |  with  a  Gaussian  density 
which  is  unimodal  by  definition.  If  the  unique  mode  of  the  density  fails  to  switch  from  the 
dominant  pre-turn  mode  to  the  dominant  post-turn  mode,  then  the  EKF  will  fail  to  hold 
track  on  the  maneuvering  target.  Much  better  performance  would  result  if  the  full  nonlinear 
problem  could  be  solved  in  real  time. 

In  this  report,  we  propose  some  new  algorithms  for  the  target  tracking  which  are  based  on 
the  spectral  nonlinear  FPS  and  the  method  of  operator  splitting.  It  has  been  demonstrated 
that  these  techniques  decrease  prediction  error  and  provide  a  more  accurate  representation 
of  the  target  bearings  as  compared  to  EKF.  On  the  other  hand,  these  algorithms  are  also 
fast,  in  that  their  calculations  need  only  0{N)  flops  per  time  step  where  N  is  the  number 
of  points  in  the  spatial  domain.  That  is,  our  approximation  of  the  optimal  nonlinear  filter 
has  an  optimal  computational  complexity  for  arbitrary  nonlinear  systems. 

In  Section  2,  the  problem  of  target  tracking  is  formulated  in  a  proper  framework  of 
nonlinear  filtering  by  using  the  basic  models  of  manuvering  targets.  Section  3  is  devoted  to 
the  development  of  the  fast  nonlinear  filters,  and  a  numerical  example  is  given  in  Section  4. 
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2  Target  Tracking  via  Nonlinear  Filtering 

Changes  in  target  orientation  is  one  of  the  most  important  sources  of  variability  of  its  sig¬ 
nature.  The  accuracy  and  the  speed  of  target  identification  algorithms  depend  crucially  on 
the  availability  and  quality  of  target  orientation  data.  In  this  Section  we  will  introduce  a 
new  model  for  tracking  dynamically  changing  orientation  via  optimal  filtering.  We  will  in¬ 
vestigate  the  flight  dynamics  equations,  develop  the  optimal  target  tracking  and  orientation 
estimation  model,  and  conduct  a  preliminary  study  of  the  angle-only  tracking  problem. 

2.1  Basic  Models  for  Maneuvering  Targets 

Let  V  —  (ui,U2,U3)  and  q  =  (91,92,93)  be  the  velocity  and  angular  velocity  of  the  target, 
respectively.  Let  /  =  (/i,/2,/3)  be  the  relative  time-derivative  of  the  velocity,  i.e.,  the  rate 
of  change  of  the  velocity  referred  to  the  body-frame  coordinates  (with  the  rate  of  change 
being  resolved  back  into  the  inertial  coordinates). 

In  this  case,  we  obtain  the  following  form  of  Euler’s  equations  (see  [24]) 

-f-  q^V2  -  92U3  =  /i, 

f’2 +  9i^’3 -93^’i  = /2,  (1) 

h  +  92^^i  -  9i^^2  =  h- 

An  important  type  of  basic  motion  in  target  maneuvering  is  constant  turn.  In  a  constant 
turn,  the  velocity  as  viewed  from  the  body-fixed  frame  is  constant.  So  /  =  0.  Also  9  is 
independent  of  time  t.  Then  equation  (1)  becomes 

vi{t)  1  To  -93  92  1  r  ui(i)  ' 

^2(f)  =  93  0  -91  V2{t)  .  (2) 

.  ^3(0  J  L  “92  9l  0  J  [  U3(f)  _ 

In  the  notation  of  cross  product  of  vectors,  this  may  be  simply  written  as  v{t)  =  9  x  v{t). 
Differentiating  both  sides  with  respect  to  t  yields 

v{t)  =  9  X  v[t) 

=  qx{qxv{t)) 

=  q{q^v{t))  -  v{t){q^q)  . 

Since  the  velocity  and  the  angular  velocity  are  perpendicular  to  each  other,  q^v(t)  =  0. 
Thus  we  end  up  with  another  constant  turn  model 

v{t)  =  -u}\(t)  ,  (3) 

where  w  =  II9II2  =  yq^q  is  the  turn  rate  or  turning  speed.  From  v  =  q  x  v  and  9  ±  u,  it 
follows  that  II9II2  =  llf’lb/lblb-  This  is  a  general  formula  for  oj  in  the  case  /  =  0. 

Remark  1.  Sometimes  constant  turn  models  are  also  called  constant  speed  turning 
models.  In  fact,  since  the  rotation  matrix  R  is  orthogonal  and  velocity  V(=  Rv)  in  the 
body-frame  coordinates  is  constant,  we  have  ||u||2  =  ||i2^y||2  =  II^IU,  i-e.,  the  speed  is 
constant. 


Remark  2.  Comparing  the  two  constant  turn  models,  we  see  that  (3)  is  formally  more 
general  than  (2);  In  (2)  both  the  direction  and  magnitude  of  the  angular  velocity  are  fixed, 
whereas  in  (3)  only  the  turning  speed  is  constant  but  the  turning  direction  may  still  change. 
Also,  model  (3)  is  decoupled  -  each  component  has  a  separate  equation  independent  of  the 
other  components,  while  model  (2)  is  not.  On  the  other  hand,  (2)  is  linear  in  q,  while  (3)  is 
nonlinear  in  uj. 

For  practical  applications,  the  above  deterministic  models  are  restrictive.  By  considering 
some  additional  white-noise  acceleration  besides  the  constant  turn,  we  obtain  from  (2)  the 
following  three-dimensional  nearly  constant  turn  model: 

Vi{t)  1  r  0  -q3  92  1  r  Vi{t)  1  r  criwi{t)  ■ 

V2{t)  =  qs  0  -qi  V2{t)  +  cr2W2{t)  ,  (4) 

.  ^3(0  J  L  “92  9i  0  J  L  ^3(1)  J  [  crsWsit)  _ 

where  cti,  <72,  and  (73  are  scaling  parameters,  and  {wi,W2,  W3)  is  a  standard  three-dimensional 
Wiener  process.  This  is  a  special  case  of  the  general  equation  (1),  where  the  angular  velocity 
q  is  constant  and  the  “relative”  acceleration  /  is  assumed  to  be  white  noise. 

Let  A  =  tk+i  —  tk  be  the  sampling  period.  Solving  the  deterministic  part  of  (4)  for  v  on 
the  interval  [4,4+i],  we  obtain  its  discretized  state  equations: 

^^i(4+i)  njCi  +  C  nin2Ci-n3S  nsniCi -f  7225  I  F  Ui(tfc)  * 

^^2(4+1)  =  nin2C'i -I- 7235  njCi  +  C  7227236'!  -  7225  U2(4) 

.  773(4+1)  J  [  723221(71  -  7225'  722723(7i  -f  72i5'  ?2^(7i  +C  J  [  773(4)  . 

ai/\wi^k+i  S  sina7A 

-f  C2/S.W2,k+i  ,  with  C  =  cos  a;  A  .  (5) 

_  (73  Aw3,fc+1  J  L  J  [  1  -  cos  07  A  _ 

Here  ni  =  9i/o7,  AtHi^it+i  =  u;,(4+i)  —  Wi{tk)  {i  =  1,2,3),  and  the  stochastic  noise  term 
has  been  simplified.  (An  exact  solution  may  be  obtained  but  will  be  much  more  complex. 
For  general  results  on  exact  discretization  of  linear  stochastic  systems,  see  [21].)  The  above 
formula  is  useful  for  describing  the  trajectory  of  a  target  when  q  =  07(721,722,723)^  is  known. 
The  deterministic  part  is  also  used  in  computer  graphics  for  rotation  models. 

Now  for  the  second  constant  turn  model,  if  we  add  some  white  noise  to  system  (3)  and 
augment  it  by  the  identities  Xi[t)  =  Vi{i)  and  Vi{t)  =  a,(t),  we  get 

Xi{t)  0  10  Xi{t)  0 

Vi{t)  =  001  Vi{t)  -b  0  criWi{t),  2  =  1,2, 3,  (6) 

hilt)  _  _  0  —07^  0  a,(t)  1 

where  x{t)  =  (a:i(t),  a;2(t),  a;3(t))  is  the  position  of  the  target  and  a  =  (01,02,(23)^  is  its 
total  acceleration  in  the  inertial  coordinates.  Solving  for  (x,-,  Uj,  a^)^  on  interval  [4,4+i],  we 
obtain  the  following  discrete-time  nearly  constant  speed  turning  model  (see  [4]): 

■  1  r  1  ^  1  r  ^i{h)  1  r  • 

Vi(tk+i)  =  0  coswA  Vi{tk)  +  5A  cTi^Wi^k+i  ,  (7) 

a, (4+1)  J  L  0  — u7sin(.(;A  coswA  J  [  0^(4)  J  .  1  . 
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where  the  noise  term  is  again  only  an  approximation  of  the  corresponding  (more  complicated) 
stochastic  integral. 

Remark  3.  The  nearly  constant  velocity  model,  the  nearly  constant  acceleration  model, 
and  the  nearly  coordinated  turn  model  (in  2D)  are  special  cases  of  the  above  two  nearly 
constant  turn  models.  More  precisely,  the  nearly  constant  velocity  model  (or  white-noise 
acceleration  model)  is  a  special  case  of  (4)  with  qi  =  q2  =  Qs  =  0,  and  the  nearly  constant 
acceleration  model  (or  Wiener-process  acceleration  model)  is  a  special  case  of  (6)  with  tu  =  0. 
The  so-called  coordinated  turn  model  in  2D  is  a  special  case  of  (2),  in  which  qi  =  q2  =  0  and 
qs  =  u: 

Vi{t)  _  0  —u)  vi(i) 

V2lt)  W)  0  ^2(0 

And,  as  in  model  (4),  if  we  a.ssume  a  white  noise  “relative”  acceleration  and  also  include  the 
equations  for  the  position,  then  we  obtain  the  nearly  coordinated  turn  model  or  essentially 
constant  speed  model  (see  [3]  [26]): 

ii(t)  '  '001  0 

X2{1)  _  0  0  0  1 

vi{t)  0  0  0  —LV 

^2(0  .  .  0  0  a;  0 

Its  discrete  time  state  equation  is  obtained  as 
■  a:i(4+i)  1  r  1  0  ^  1  r  ]  r  iA  0  1 

X2(4+i)  ^  0  1  ^  X2{h)  ^  0  [  criAu;i,,+i  1 

i’i(4+i)  0  0  coscuA  —  sincuA  Vi{tk)  1  0  <72Atr;2,fc+i  ’  ^ 

_  ■^2(4+1)  J  L  0  0  sintuA  coswA  J  [  ^2(4)  J  L  0  1  J 

2.2  Estimation  of  Target  Orientation 

The  purpose  is  to  estimate  the  orientation  angles  of  the  target.  This  is  achieved  by  using 
the  basic  models  derived  above  and  the  Euler’s  relation  between  the  angular  velocity  and 
the  orientation  angles. 

One  approach  using  the  first  nearly  constant  turn  model  has  been  developed  in  earlier 
work  of  this  project  (and  in  particular  in  Phase-I  work  of  this  project).  Here  we  present 
another  approach  using  the  second  nearly  constant  turn  model. 

To  estimate  the  turn  rate,  we  assume  cu  is  a  diffusion  process  satisfying 

uj{t)  =0(^-1-  7i^ci;(t)  +  a^^w^{t). 

Then,  from  this  and  (6),  we  obtain  a  3-D  system 

v{t)  =  a{t), 

a{t)  =  -u‘^{t)v{t)  -h  CTaWait),  (10) 

Uj{t)  =  a^-\-  7i^u;(f)  -f-  a^w,^{t), 

where  a^,  and  cTa  are  constants,  and  Wa{t)  and  w^{t)  are  independent  standard  Wiener 

processes. 


'  a:i(t)  ■  '00' 

^2(0  I  0  0  (XiWi{t) 

Vl{t)  1  0  Cr2W2{t) 

.  L  ^^2(4  J  L  0  1  J 
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As  for  the  measurements,  aissume  the  position  x{t)  of  the  target  is  observable  at  discrete 
time  t  =  tk-  Then  x(4+i)  —  x{tk)  is  also  observable.  From  (7),  we  have  the  following 
observation  equation: 

sin(a;(4)A)  +  J.  x  a(4)(l  -  cos(a;(4) A))  +  e,v^{k)  ,  (11) 

where  {n2(fc)}  is  a  sequence  of  independent  Gaussian  random  variables  of  zero  mean  and 
unit  variance,  and  Cz  is  the  standard  deviation  of  the  noises. 

Now  we  have  established  model  (lO)-(ll)  for  the  velocity,  acceleration  and  turn  rate  of 
the  target.  To  obtain  a  satisfactory  estimate,  we  need  a  good  filtering  algorithm.  This  will 
be  discussed  in  Section  3. 

From  the  velocity,  acceleration,  and  turn  rate,  the  angular  velocity  can  be  calculated. 
After  obtaining  the  angular  velocity,  we  can  calculate  the  target  orientation  as  follows.  Let 
0(f)  =  (01  (f),  02(^),  be  the  Eulerian  orientation  angles  of  the  target  at  time  t  (resolved 
in  the  fixed  coordinate  system).  Then  the  components  (^i,  ^2)  93)  of  the  angular  velocity  can 
be  expressed  in  terms  of  (0i,02,03)  (see  [24]): 

qi  =  01  cos  02  cos  03  —  02  sin  03, 

q2  =  01  cos  02  sin  03  +  02  cos  03,  (12) 

93  —  sin  02  +  03. 

Solving  for  (0i,02,^^3)  gives 

01  =  {qi  cos  03  +  (j2  sin  4>z)  sec  02, 

^2  =  -qi  sin  03  +  q2  cos  03,  ( 13) 

<i>3  =  {qi  cos  03  +  (J2  sin  03)  tan  02  +  qs- 

Thus,  to  obtain  the  conditional  expectation  and  conditional  covariance  of  the  orienta¬ 
tion  angles,  one  can  approximate  (13)  directly.  These  are  deterministic  ordinary  differential 
equations,  which  can  be  easily  solved  by  using  Runge-Kutta  method,  for  example.  (Another 
approach  to  obtain  the  expectation  and  covariance  of  the  angles  is  to  first  replace  the  dynam¬ 
ics  equations  by  an  extended  system,  including  the  angles,  and  also  rewrite  the  observation 
equation  in  terms  of  all  the  new  state  variables.  Then  consider  a  higher-dimensional  filtering 
problem  based  on  the  new  equations.  This  will  give  the  estimation  of  both  the  angular 
velocity  and  the  orientation  angles  at  the  same  time.) 

2.3  Target  Tracking  with  Angle-only  Measurements 

There  are  military  situations  in  which  it  is  desirable  to  estimate  the  position,  velocity  and 
perhaps  acceleration  of  a  target  from  measurements  of  angle  but  not  range.  A  well-known 
example  is  the  determination  by  a  submarine  of  planar  position  and  velocity  of  a  ship  from 
passive  sonar  measurements,  because  the  submarine  commander  does  not  want  to  reveal  his 
presence  by  pinging.  In  air  warfare,  a  fighter  defending  against  a  raid  may  wish  to  launch 
a  missile  against  a  jammer  at  unknown  range,  but  should  not  do  so  unless  the  jammer’s 
position  and  velocity  can  be  estimated.  A  more  recent  problem  is  estimation  of  target 
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position,  velocity  and  acceleration  in  three  dimensions  from  angle  measurements  only,  either 
with  a  passive  IR  receiver  or  a  jammed  radar  receiver  on  a  missile,  in  order  to  utilize  optimal 
guidance. 

There  appear  to  have  been  two  main  approaches  to  angle-only  tracking:  (a)  tracking 
based  on  the  Extended  Kalman  Filter  in  Cartesian  coordinates  (see  [1]  [9]);  (b)  reformulation 
of  the  problem  in  terms  of  modified  polar  coordinates  with  subsequent  application  of  EKF 

([14]). 

Unfortunately  neither  one  of  the  above  is  satisfactory.  The  main  limitation  on  EKF  ap¬ 
plication  to  angle-only  tracking  is  the  nontrivial  nonlinearity  of  the  latter.  Its  mathematical 
model  can  be  described  as  follows  (for  the  sake  of  simplicity  we  describe  here  the  2D  case 
and  a  3D  example  will  be  given  in  Section  4): 

Signal: 

dxi{t)  =  bi{xi{t),X2{t))dt  -I-  o-idwi{t),  a:i(0)  =  a;?, 
dx2{t)  =  b2{xi{t)^  X2{t))dt  -b  a2dw2{t),  0:2(0)  =  x”; 

Observation: 

dy{i)  =  arctan(x2(t)/xi(t))dt  -1-  crdv[t),  t/(0)  =  0, 

where  Wi{t),W2{t)  and  v{t)  are  independent  Wiener  processes. 

The  modified  polar  coordinates  (MFC)  introduced  by  Hoelzer  et  al.  [14]  are  designed  to 
reduce  the  nonlinear  observations  to  linear.  Unfortunately  in  doing  so  MFC  also  transforms 
the  signal  process.  Unless  the  latter  is  of  very  simple  nature,  the  MFC  transform  makes  the 
signal  system  practically  intractable. 

Much  more  perspective  approach  to  angle-only  tracking  is  based  on  optimal  nonlinear 
filtering.  The  optimal  nonlinear  filtering  theory  allows  to  compute  the  posterior  density 
function  p{t,x)  of  the  signal  process  x{t)  =  (xi{t),X2(t))  given  observations  y{s),  s  <  t. 
The  posterior  density  function  is  defined  by  p{t,x)  =  u{t,x)l  f  u{t,x)dx  where  the  so  called 
unnormalized  filtering  density  u{t,  x)  is  a  solution  of  the  Zakai  equation 


du(t,  x)  =  jCu(t,  x)dt  -h  h{t)dy{t),  tt(0,  x)  =  p(0,  x). 


where 


/  d{biu)  d{b2u)\ 

\  dxi  dx2  )  ' 


The  optimal  estimate  of  the  signal  (xi(t),  X2(t))  is  given  by 


Xi{t) 


f  XiU{t,Xi,X2)dxidx2 

J  u{t,Xi,X2)dXidx2  ’ 


For  a  long  period  of  time  practical  application  of  nonlinear  filtering  has  been  strained 
by  numerical  difficulties  related  to  on-line  solution  of  the  Zakai  equation.  Recently  this 
problem  was  resolved  with  the  introduction  of  a  spectral  approach  to  nonlinear  filtering  by 
B.L.  Rozovskii  and  his  co-workers  (see  [17]).  The  spectral  approach  proposed  in  these  works 
is  based  on  Wiener  Chaos  expansion.  An  important  feature  of  this  expansion  is  that  it 
separates  observation  (t/(s),  s  <  t)  and  parameters  (61, 62,1X1,  (72, cr)  and  h  =  arctan(^)  in 
Zakai ’s  equation.  The  numerical  algorithm  for  solving  the  Zakai  equation  based  on  the  spec¬ 
tral  approach  splits  into  two  parts:  “deterministic”  and  “stochastic”.  The  time  consuming 
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computation  of  the  deterministic  part  can  be  shifted  off-line.  The  stochastic  part  involving 
the  observation  process  y{t)  is  computationally  simple  and  can  be  performed  in  real  time. 

In  the  course  of  this  project  we  compared  an  angle-only  tracker  based  on  a  spectral 
algorithm  for  nonlinear  filtering  prediction  and  smoothing  with  a  standard  EKF  tracker.  We 
considered  several  practically  important  cases,  in  particular:  (1)  the  target  exhibits  special 
evasive  maneuvers;  (2)  lack  of  prior  information  about  position  and/or  speed  of  the  target. 
In  all  cases  the  performance  of  the  nonlinear  angle-only  tracker  was  superior  to  the  EKF 
tracker. 


3  Fast  Nonlinear  Filters  of  Linear  Complexity 

Our  objective  here  is  to  develop  recursive  numerical  algorithms  for  computing  the  optimal 
filter  in  which  the  on-line  computation  is  as  simple  as  possible;  in  particular,  the  number  of 
computer  operations  at  each  time  step  should  be  proportional  to  the  number  of  grid  points 
where  the  filtering  density  is  numerically  defined.  The  starting  point  in  the  derivation  is 
the  equation  for  the  unnormalized  filtering  density  in  the  general  nonlinear  model,  and  the 
approach  is  based  on  the  technique  known  as  operator  splitting. 

Since  in  most  cases  the  observational  measurements  are  only  available  at  discrete  time 
moments,  in  this  section  we  consider  dynamic  systems  with  discrete  observations. 

Two  algorithms  are  developed;  one  is  described  in  subsection  3.2  and  the  other  in  sub¬ 
section  3.3.  The  computational  complexity  of  both  algorithms  is  indeed  0{N)  where  N  is 
the  number  of  points  in  the  spatial  domain.  Their  approximation  errors  are  also  estimated. 

3.1  The  Nonlinear  Filtering  Problem 

Now  consider  the  dynamical  system  described  by  the  stochastic  differential  equation 

dXt  =  b{Xt)dt  +  a{Xt)dWt,  t>0, 

^0  ~  7ro(x), 

and  the  discrete  observations  given  by 

zk  =  h{Xt,)  ^  e{Xt,)Vk,  A;  =  l,2,---,  (15) 

where  ttq  :  jK"  ->  iR  is  the  initial  density,  b  :  FT,  and  h  :  IR^  are  known  vector¬ 
valued  functions,  a  :  FC  — >■  and  e  :  JR”  jfirnxm  matrices  with  known  function 

entries,  {Wt}t>o  is  a  standard  r-dimensional  Brownian  motion,  {14}fc>i  is  a  standard  di- 
dimensional  white  Gaussian  sequence,  and  tk  =  kA  (A  >  0).  Xq,  {Wf}  and  {14}  are 
assumed  to  be  independent,  and  functions  6,  h,  cr,  e  and  ttq  are  assumed  to  be  smooth  enough 
(satisfying  certain  regularity  conditions,  see  [19] [23]). 

Let  /  =  f{x),x  G  iR",  be  a  measurable  scalar  function  such  that  lE\f{Xt)\^  <  oo  for  all 
t  >  0.  Then  the  filtering  problem  for  (14)-(15)  can  be  stated  as  follows:  find  the  minimum 
variance  estimate  of  /(Ai^)  given  the  measurements  •  •  • ,  z^.  This  estimate  is  called  the 
optimal  filter  and  is  known  to  be 

fk  =  E[f{Xt,)  \zi,---,Zk  . 
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For  computational  purposes,  the  optimal  filter  can  be  characterized  as  follows. 

Denote  by  Tt  the  solution  operator  for  the  Fokker-Planck  equation  corresponding  to  the 
state  process;  in  other  words,  u{t,  x)  =  Tt^p{x)  is  the  solution  of  the  equation 


du{t,  x) 
dt 


\  E  ({a{x)a{xf)^,u{t,  x)) 

n  Q 

-  ^  —  {h^{x)u{t,  x)),  t>  0, 

U=l 


u(0,  a:)  =  <p(a:). 


(16) 


m 

where  (cr(x)cr(x)^)^^  =  is  the  ^-th  row  and  z/-th  column  entry  of  cr(a:)cr(x)^, 

i=l 

and  hv{x)  is  the  i/-th  component  of  b[x). 

Next,  define  the  unnormalized  filtering  density  pk{x),  for  x  G  and  A;  >  0,  by 


Po{x)  =  TToix), 

Pk{x)  =  ak{x)TAPk-iix), 

where 

afc(a;)  =  exp  -  h{x)f{€{x)e{xf)~\zk  -  /i(a;))| ,  A:  =  1, 2,  •  •  • . 

Then  the  optimal  filter  fk  can  be  written  as  follows  [15]: 


(17) 


i  ^  Pk{x)f{x)dx 
^  fR’^Pk{x)dx 


(18) 


3.2  Splitting  of  Convection  and  Diffusion  Terms 

To  compute  the  unnormalized  filtering  density,  a  fast  Fokker-Planck  solver  is  needed.  In 
previous  work  related  to  this  project,  two  methods  have  been  developed:  a  method  based 
on  the  spectral  separation  approach  ([18]),  and  a  method  based  on  the  finite  element  ap¬ 
proximation  ([19]  [23]).  In  this  subsection  we  present  another  method  which  is  based  on  the 
operator-splitting  technique  ([16][20]).  Similar  ideas  have  been  used  in  [13]  for  solving  the 
Zakai  equation  arising  from  image  filtering. 

We  first  assume  that  in  the  noise  term  of  (14),  n  =  r  and  the  covariance  matrix  a  is 
diagonal  and  constant:  crf^^(x)  =  5^,ya^.  Then  the  Fokker-Planck  equation  (16)  becomes 


du{t^x)  1  " 

dt  2  E  dxl 

u(0,  x)  =  <p(x)  . 


{alu{t,x)) -'^—{h^{x)u{t,x)y  t>0. 


U=\ 


dx 


Its  solution  can  be  expressed  as 

Tt^{x)  =  E[ip{^^)exp  {  -  (V  •  6)(^J)ds}]  , 


(19) 
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where  (f  is  a  stochastic  process  satisfying 


—  —b{^f)dt  +  di&g{a^)dWt  , 
=  x]  =  1  . 


To  proceed  the  splitting  of  convection  and  diffusion  terms,  denote  by  T/  and  Tf  the 
solution  operators  of  the  equations 

du{t,x)  A  ^  /'i  /  ^  /.  nA  .  n 

=  -  E  ^  (^)^(a:)u(t,  x) j ,  t>  0, 

w(0,x)  =  (f{x)  , 


and 


du{t,x)  1  ^  /  2  /-j  \\  j  n 


u(0,x)  =(p{x)  , 

respectively.  In  fact,  the  two  solution  operators  can  be  expressed  explicitly  as 

rt 


and 


r<V(x)  =  v?(?7r)exp  {  -  ^  (V  •  6)«)ds}  , 

(7rt)-"/2  r  .  "  (x,-  -  y,)2 


T,V(x)  =  Ei^iCn]  =  /  3  exp  {  - 

ai-'-Un  Jr^  ^ 

where  processes  T]f  (deterministic)  and  Q  satisfy 

dVit  =  Vo=x  , 


2alt 


}<p{y)dy , 


(20) 


(21) 


and 

dCt  =  diSig{a^)dWt,  Co  =  ^  • 

Then  it  can  be  shown  that  the  following  approximation  formulas  hold: 

TaV  =  TiTlv  +  0(A) ,  (22) 

TaV  =  TlTiTlf  +  0(d})  .  (23) 

2  2 

Therefore,  instead  of  solving  the  original  Fokker-Planck  equation,  we  only  need  to  com¬ 
pute  (20)  and  (21).  To  compute  (20),  we  only  need  to  solve  an  ordinary  differential  equation. 
To  compute  (21),  we  only  need  to  integrate  over  a  small  area  near  point  x,  especially  when 
the  noises  are  small.  In  the  case  of  large  noises,  the  multi-grid  method  ([12])  can  be  used  to 
achieve  linear  computational  complexity. 

Remark  1 .  The  Trotter’s  semigroup  approximation  theorem  guarantees  the  convergence 
of  the  above  approximations  (22)  and  (23)  (but  it  does  not  provide  any  error  estimate).  In 
our  notation  it  implies 

T,v  =  ^lim  (r,%r,%)  V . 
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Remark  2.  To  see  the  difference  between  the  exact  solution  and  the  approximate  solu¬ 
tion  (by  splitting),  we  note  that 


r/r,v(x)  =  iE[T,v(cf) 


=  E[ip  (?7t(C<(a:)))  exp  {  -  ^  ( ^  •  b){Tjs{Ct{x)))ds}] , 


T^TfTM^)  =  r/rf^(77|)exp{  -  •  6)(r7f)ds} 

=  E[ip  (?7i(G(»7|W)))  exp  I  -  Jj{W  ■  6)(7/,(0(?/i(a;))))ds}j  exp{-  Jj{V  ■  b){r]s{x))ds] 
=  e[(p  {r)i{Ctivi{x))))  exp  {-  b)r]s{x))ds  - 


where  for  obvious  reasons  we  have  written  r]f  =  r]t{x)^  (f  =  Ct(x),  and  rjA  =  (^t{x,t'),  the 
last  being  the  process  rjt  starting  from  the  point  x  at  time  t'.  Comparing  these  results  with 
(19),  we  find  that  in  effect  our  approximation  is  to  split  the  process  into  two  simpler 
processes:  a  deterministic  process  rjt  and  (up  to  a  constant)  a  Wiener  process 

In  general,  if  the  covariance  matrix  a  —  (r[x)  is  not  constant  (but  still  assumed  to  be 
diagonal  to  simplify  expressions,  i.e.,  cr^,y(a;)  =  S^^a^lx)),  then  we  rewrite  the  Fokker-Planck 
equation  in  the  following  form: 


du{t,x)  1  ^ 

2  dx 


dt 


U- 

n 


+  E 


u-\ 


dxi, 

(^(a^{x)- — a^{x)  -  b^{xfju{t,x)y 


and  then  split  it  into  the  following  two  equations: 


du{i,  x) 
dt 


i/=i 


by(x)ju(t,x)), 


and 


du(t,  x) 
dt 


both  of  which  can  be  solved  in  linear  computational  complexity. 

In  this  general  situation,  the  three  solution  operators  T),  T/,  and  Tf  can  also  be  expressed 
via  certain  stochastic  (or  even  deterministic)  processes.  Indeed,  similar  to  formulas  (19),  (20), 
and  (21),  we  have 


Tt^{x)  =  E  [97(^f )  exp  {  ^  V  •  (a  -  fe)(Ods}] , 
T^^{x)  =  v3(»7f)  exp  {  ^  V  •  (d  -  6)(?7j)ds}, 
r,V(x)  =  E[<p{C,t)]  , 
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where  processes  ,  r/f,  and  Cf  satisfy 


dff  =  {2a(S)  -  mm  +  d\^(a,(S))dW„  fj  =  T, 

drif  =  -  b{ttfm<  nS  =  x, 

and 

<<cr  =  +  <liag(a„(cr))<itr.,  Co'  =  i:, 

respectively.  Here  we  have  denoted  by  d(y)  the  vector  with  components  a^(j/)  ^0^,(1/)  (1/  = 

1 ? ' '  ’  ?  • 

Remark  3.  The  algorithm  by  splitting  the  convection  (or  drift)  term  from  the  (pure) 
diffusion  term  as  discussed  above  has  a  further  advantage  that  much  of  the  computation  in 
(20)  and  (21)  or  their  generalizations  can  be  performed  before  the  observations  are  available. 
This  pre-calculation  can  save  the  on-line  cost  and  further  speed  up  the  real  time  performance. 


3.3  New  Alternating  Direction  Implicit  Schemes 

In  the  previous  subsection  we  used  the  operator-splitting  technique  to  split  the  whole  process 
into  a  convection  process  and  a  diffusion  process  (without  “drift”).  In  this  subsection  we 
will  apply  the  operator-splitting  technique  to  split  a  multi-dimensional  process  into  several 
one-dimensional  processes,  since  the  filtering  densities  for  1-D  processes  can  be  computed 
easily  (with  linear  computational  complexity).  And  this  is  the  idea  of  alternating  direction 
implicit  (ADI)  schemes  ([11]  [16]  [20]  [27]). 

Here  we  present  two  new  splitting  schemes  based  on  the  discretization  method  de¬ 
scribed  in  [25];  one  is  similar  to  the  Dyakonov  scheme  ([7]  [8])  and  the  other  is  similar  to 
the  Peaceman-Rachford  scheme  ([22]).  We  first  discuss  the  2D  case  in  details  and  then 
generalize  the  results  to  higher  dimensions.  Our  higher  dimensional  generalization  of  the 
2D  Peaceman-Rachford-like  scheme  is  much  simpler  than  the  usual  ADI  modifications  of 
Peaceman-Rachford  scheme  in  higher  dimensions  ([5]  [6]  [16]  [20]). 

First,  let  us  consider  the  two-dimensional  convection-diffusion  equation 

Ut  =  auxx  +  huyy  -b  cux  +  duy  +  eu,  in  (0,  T]  x  0,  (24) 

and  the  initial-boundary  value  conditions 

u(0,  X,  y)  =  u°(x,  y),  (x,  y)  e  D,  ,  . 

u(f,x,y)  =  0,  t  G  (0,r],(x,y)  e  5D,  ''  ’ 

where  0  =  (a,/?)  x  (7, 5), a  <  (3,^  <  6,T  >  0]  a,b,c,d,  and  e  are  known  smooth  functions 
with  bounded  derivatives,  defined  in  Dx  :=  [0?^]  x  0;  G  (7(fi)  satisfies  the  continuity 
condition  u°(x,y)  =  0,V(x,y)  G  dCl.  For  simplicity,  we  assume  min(o,6)  >  A  >  0  (A  a 
constant),  and  e  <  0. 

Let  X,  :=  a-bfAx  (0  <  i  <  Nx),  yj  :=  q-bjAy  (0  <  j  <  Ny),  and  4  :=  kAt  {0  <  k  <  M) 
be  a  partition  of  fix?  with  Ax  :=  Ay  :=  and  At  :=  For  any  v  G  C'(fix))  denote 
v’^-  :=  t;(4,  Xj,  yj)  {0  <  i  <  Nx,  0  <  j  <  Ny,0  <  k  <  M). 
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We  approximate  problem  (24)-(25)  by  the  following  difference  equations: 


At 


^k+^D+yD-y\  ( 

^3.2  ij  Ay2  y  V  2  J 

‘'■-+  V  2Ax  j  2Aa:  ) 

,  jk+l  f  D+yV^j^^  +  D-yV^j\  ,  f  D-yV^t^  +  D+yvfj\ 

+  ^  2Ay  ;  +  ^  2Ay  ) 

S1±A\ 

2  r 


+  e 


^+1  r 


f  =  !,•••,  fV^-1,  i  =  l,---,Arj,-l,  fc  =  0,1,- ••,M-1; 


=  U°( 

=  0,1,-" 

'  )  A^a:)  j  —  0)  I)  ■  ■  ■ 

1  ^y> 

=  0, 

<,i  =  0. 

i  =  1,- 

Ay-1,  k=l, 

=  0, 

0 

II 

f  =  1,  • 

,Nx-l,  k  =  1, 

(26) 


(27) 


where  D±x  and  D±y  denote  the  forward  and  backward  difference  operators  in  the  x  and  y 
directions,  respectively,  are  defined  as 


%•.+  ■ 


k-j-  h 

r.  .  ^ 


1  /c+1 


=  max(0,cf/^)  = 

L-.,  >0]  V  5  o  /  2 

*+211  ■  /n  *+2\ 

=  . . . 


Mk 


and  similarly  for  ^  •  It  can  be  shown  (see  [25])  that  the  truncation  error  of  scheme  (26) 
is  0(^{At  +  Ax  +  Ay)^^  and  that  the  scheme  is  unconditionally  stable. 

Multiplying  both  sides  of  (26)  by  At,  we  can  rewrite  system  (26)-(27)  into  the  following 
matrix-vector  form: 


( 


I - ^(^S.l  +  ^y,l)^  4-  — (Aj;,0  -f  Ay^o)  ) 


(28) 


where  matrices  Axfi,Ax^i,  Ayfi,Ay^i  axe  “tridiagonal”:  all  but  two  of  their  off-diagonals  are 

zero  and  the  two  nonzero  off-diagonals  are  of  the  same  distance  from  the  diagonal;  and 

jt+i 

each  of  these  matrices  contains  e,j  ^/2  in  the  diagonal  entries.  (In  general,  these  coefficient 
matrices  depend  on  the  time  step  k,  but  we  have  ommitted  the  superscript  A;  4- 1  to  simplify 
the  discussions  that  follow.) 

As  usual,  equation  (28)  is  numerically  much  more  difficult  than  its  one-dimensional  ana¬ 
logue  (see  [25]).  As  a  resolution  to  overcome  this  difficulty,  we  replace  (28)  by  the  following 
“one-dimensionalized”  approximation 


(29) 


Now  we  only  need  to  solve  a  sequence  of  tridiagonal  systems  at  each  time  step. 
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Since  scheme  (29)  differs  from  scheme  (28)  only  by  the  term 


which  is  easily  checked  to  be  of  order  0{At^{At  +  Ax  +  Ay))  (meaning  a  difference  of 
truncation  error  by  0{At^  +  At  Ax  +  At  Ay)),  the  accuracy  of  scheme  (29)  remains  the  same 
order,  that  is,  0{{At  +  Ax  +  At/)^). 

An  alternative  approach  is  to  use  the  following  alternating  direction  implicit  (ADI) 
scheme: 


(/+ Y-4,,o)t.‘, 

=  (/  +  ^A.,o)  . 


(30) 


To  obtain  the  truncation  error  of  this  Peaceman-Rachford-like  approximation,  we  notice 
that  equations  (29)  and  (30)  can  be  written  as  follows: 


When  compared  with  (31),  the  term  (/  — is  replaced  by  {I+^Axfi){I— 
in  (32).  But  both  of  these  are  solving  the  same  “one-dimensinal”  convection- 
diffusion  equation  with  the  same  order  of  approximation  error.  So  the  accuracy  of  scheme 
(32)  is  the  same  as  that  of  (31),  i.e.,  0{{At  +  Ax  +  AyY). 

Note  that  if  traditional  Crank-Nicholson  discretization  is  used,  i.e.,  if  Ax,].  =  Ax,o  and 
Ay,i  =  Ay,o,  then  (29)  and  (30)  become  Dyakonov  and  Peaceman-Rachford  schemes,  respec¬ 
tively.  In  this  case,  the  two  schemes  are  in  fact  equivalent;  i.e.,  in  two  dimensions,  Dyakonov 
scheme  and  Peaceman-Rachford  scheme  are  just  two  different  forms  of  the  same  approxima¬ 
tion.  On  the  other  hand,  with  our  discretization,  schemes  (29)  and  (30)  generally  are  not 
equivalent. 

The  above  two-dimensional  schemes  can  be  generalized  to  higher  dimensions.  Indeed, 
similar  to  (31)  and  (32),  we  may  use  one  of  the  following  two  splitting  schemes  for  three- 
dimensional  homogeneous  problems: 


(/  +  y.4.,„)  (/  +  ^A,,„)  (/  +  ^A,,„)  /  i 

And  in  the  four-dimensional  case,  we  use  one  of  the  following  schemes: 


(33) 


(34) 
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(35) 


The  same  arguments  as  in  the  2D  case  can  be  used  to  show  that  all  these  schemes  are  second 
order  accurate.  And  it  can  also  be  shown  that  all  these  schemes  are  absolutely  stable. 

Finally,  we  make  two  remarks  about  the  above  ADI  schemes. 

1.  Whereas  the  order  of  the  alternating  directions  or  variables  in  the  Dyakonov-like 
schemes  (31),  (33)  and  (35)  does  not  affect  the  second-order  accuracy  (up  to  higher  order 
terms  as  long  as  the  explicit  and  implicit  schemes  are  applied  separately),  an  impropriate 
order  of  alternating  directions  in  the  Peaceman-Rachford-like  schemes  (32),  (34)  and  (36) 
would  lead  to  first-order  accuracy.  For  example,  in  the  3D  case,  in  (33)  we  may  also  use  the 
order  y,z,x,z,y,x  instead  of  x,y,z,  z,y,x.  But  in  (34),  we  can  only  change  the  order  to 
something  like  y,  z,  x,  x,  z,  y;  the  second  half  must  be  in  the  inverse  order  of  the  first  half. 

2.  Our  higher  dimensional  generalizations  (34)  and  (36)  are  simpler  than  the  usual  higher 
dimensional  ADI  modifications  of  Peaceman-Rachford  scheme,  even  for  Crank-Nicholson 
discretization  with  Ax^fi  —  Ax^,i  —  Ax^i'iv- 


4  Numerical  Example 


To  illustrate  the  performance  of  the  above  algorithm,  let  us  consider  the  following  three- 
dimensional  tracking  problem  with  angle-only  observations 


Signal: 


Observation: 


Mt)  =  bnXiit)  -h  b^2X^{t)  +  <7iiWi(0, 
X2{t)  =  b2  +  Cr22W2{t), 


Xsit)  =  63-1-  <T33W3(t), 

z,(k)  =  sign(*(( J)  arccos  +  Mk), 

Z2(k)  =  arcsin 


where  611  =  —200,  612  =  50,  62  =  —1/2,  62  =  —1/4,  cth  =  0.045,  <722  =  0.023,  <733  =  0.012, 
tk  =  kA,  A  =  0.01,  Vi{k)  ~  N{0,  .64^),  V2{k)  ~  A^(0,  .36^),  and  the  initial  target  state 
(A’i(O), ^2(0), ^3(0))  has  joint  density  7ro(a:i,a:2,a:3)  =  A exp(— 100(a:i  —  0.45)^  —  121(a:2  — 
0.42)^  —  64(0:3  —  0.20)^),  Co  being  the  normalizing  constant. 

In  our  experiments,  we  took  T  =  2.0  (so  M  =  200),  S  =  [—0.95,0.75]  x  [—0.60,0.60]  x 
[—0.30,0.30],  x\  =  —0.95  -|-  iAxi  {0  <  i  <  rii),  x{  =  —0.60  jAx2  {0  <  j  <  702),  0:3  = 
—0.30  4-  A:Aa;3  (0  <  /?  <  ns),  Axi  =  Ax2  —  Axs  =  and  [ni,n2,n3]  =  [40,30,20] 
or  [80,60,30]. 
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Results  of  typical  tracking  steps  are  shown  in  the  figures  at  the  end  of  this  report,  where 
the  true  initial  target  position  was  (Xi(0),  ^2(0)5 -^3(0))  =  (0.45,0.5,0.25).  Programming 
source  code  for  this  and  other  examples  is  also  attached. 

In  conclusion,  the  filtering  algorithms  proposed  in  this  report  are  based  on  the  exact 
optimal  filter  and  so  apply  to  systems  with  high  nonlinearity  and  different  noise  levels.  On 
the  other  hand,  the  new  algorithms  are  fast,  in  that  their  calculations  have  linear  comlexity 
per  time  step. 
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Dynamics  without  noise 


-0.2 


Observations  without  any  noise  in  the  system 


Angle-Only  Observations 


Marginal  density  for  X1  (0) 


Marginal  density  for  (X1(0),X2(0)) 


Marginal  density  for  X1  (1 .0) 


Marginal  density  for  (X1  (1 .0),X2(1 .0)) 


Marginal  density  for  X1  (1 .5) 


Marginal  density  for  (X1  (1 .5),X2(1 .5)) 


Marginal  density  for  X1  (2.0) 


Marginal  density  for  (X1  (2.0),X2(2.0)) 


M.-L.  estimate  for  (X1  ,X2,X3)  Density  contour  for  (X1  (2.0),X2{2.0)) 


C.Rao 

run.m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  New  operator  splitting  method  for  solving 

%%  dX  =  b(X)dt  +  sigmaa(X)dW,  X(0)  =  mO 

%%  z(k)  =  h{X(t_k))  +  delta(k)v(k) 

%%  Chuanxia  Rao 

%%  January  1996 

%%  January  1997 

%%  May  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  This  is  the  main  program. 

%% - 

%%  Dependencies: 

%%  initialize. m 

%%  moreinit.m 

%%  generate. m 

%%  online,  m 

%%  plotting. m 

%% - 

%%  Output: 

%%  graphs 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


clear; 


count  =  flops; 
cpu  =  cputime; 

initialize;  %%  initial  data 

moreinit;  %%  further  initialization 

generate;  %%  to  generate  processes  X  and  Z 

flops_off  =  flops  ~  count 
cpu_off  =  cputime  -  cpu 


figure; 

xyz  =  X ( 1 , : ) ; 
plotting (xyz ,  pk) ; 
kstart=l ;  kstop=Kmax; 

%kstart=Kmax+l ;  kstop=Kf inal-1 ; 

count  =  flops; 
cpu  =  cputime; 

%%global  time ; 

for  ktime  =  ks tart : ks top, 

kl  =  ktime  +  1; 

%%  time  =  to  +  dt*ktime; 

fprintfd,  '%c' ,  '  step  '); 

fprintfd,  '%d\n',  ktime); 

online; 

kplot  =  ktime/plot_step; 
if  ktime==l  |  kplot  ==  fix{kplot), 
xyz  =  X(kl, : ) ; 
plotting (xyz ,  pk) ; 

end 

end 

count  =  flops  -  count; 

cpu  =  cputime  ~  cpu; 

cpu_onl  =  cpu  /  (kstop-kstart+1) 

flops_onl  =  count  /  (kstop-kstart+1) 


C.Rao 

initialize.!!! 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  Initial  data  for  the  cont.-discr.  model: 

%%  dX  =  b(X)dt  +  sigmaa*dW 

%%  X(0)  -  N(mO,varO) 

%%  z(k)  =  h(X{t_k))  +  delta*v(k) 

%% - 

%%  Chuanxia  Rao 

%%  Februay  1996 

%%  Februry  1997 

%%  June  1997 

%% - 

%%  Called  in: 

%%  all  other  parts 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

global  dim  dim_obs 
global  nx  ny  nz 
global  xa  ya  za 
global  xb  yb  zb 
global  mO  varO_inv 

global  substeps 
global  small  %large 
global  resol_vert 
global  method 

global  denom_obs 
global  denom 
global  XX  yy  zz 
%global  nxpl  nypl  nzpl 


disp('  Initialization 
dim  =  3 ; 
dim_obs  =  2 ; 

%%  Part  1. 

Tfinal  =  2.0; 

Kfinal  =  256;  Kt  =  Kfinal; 
Kmax  =  128;  %64;  %256; 

%Kmax  =  Kfinal  -  1; 

%Kmax  =  Kfinal  /  4; 


nx  = 

72; 

%60; 

%40; 

%80 

ny  = 

60; 

%45; 

%30; 

%60 

nz  = 

30; 

%24; 

%20; 

%30 

xa  = 

-0. 

95; 

xb  = 

0 

.75 

ya  = 

-0. 

60; 

yb  = 

0 

.60 

za  = 

-0. 

30; 

zb  = 

0 

.30 

%xa 

=  [- 

0.95, 

-0.6 

-0 

%xb 

=  [0 

.75, 

0.6, 

0 

.3] 

%Nx 

=  [nx,ny. 

nz]  ; 

%Lx 

=  [3 

,2,1] 

7 

Lx  = 

[1, 

1,1]  : 

%%  bandwidth  of  local  propogation: 

%%  It  depends  on  (dx*dx)  /  (dt*sigmaa"^2 )  . 


plot_step  =  2;  %1; 
resol_vert  =  0.01; 


%%  plot  in  every  plot_step  steps 
%%  for  plotting  target  position 


%%  threshold  for  positive  probability 


substeps  =  1; 
small  =  le-16; 

%large  =  -log (small); 


C.Rao 
initialize.m 

%%  Part  2 . 

%sigmaa  =  [0 . 0,  0 . 0,  0 , 0] ;  %%  coefficients  in  signal  noise 

sigmaa  =  [0,045,  0.023,  0.012]; 

%delta  =  [0.0,  0.0];  %%  coefficients  in  observ.  noise 

%delta  =  [0.24,  0.12] ; 
delta  =  [0.64,  0.36]; 

mO  =  [0.45,  0.42,  0.20];  %%  initial  mean 

var0_inv  =  [100,  121,  64] ;  %%  initial  variance,  its  inverse 

%var0_inv  =  [82,  100,  64]; 

%X0  =  [0.75,  -0.5,  -0.25] ; 

XO  =  [3.674/2,  1./2,  1./4] ; 


zo 

=  [0,0]; 

%% 

Part  3 . 

method  =  11; 

%01: 

cdsplit 

%% 

convection-diffusion 

splitting 

(cd) 

%02: 

desplit 

%% 

convection-diffusion 

splitting 

(dc) 

%11; 

edesplit 

%% 

convection-diffusion 

splitting 

(ede) 

%12: 

dedsplit 

%% 

convection-diffusion 

splitting 

(ded) 

%21: 

centrall 

%% 

central  difference  + 

backward  Euler 

%22: 

central2 

%% 

central  difference  + 

C-N  (trapezoid) 

%31: 

upwind 

%% 

1st  order  upwind  +  backward  Euler 

%32: 

new2nd 

%% 

new  "upwind"  (2nd  order  in  space  &  time 

%%  and  resulting  in  d.d.  tridiagonals) 


%% 

%%  Further  initialization; 
%%  —  no  changes  needed  -- 

%%  , 

%%  Part  4. 
dt  =  Tfinal/Kfinal; 
dx  =  (xb-xa) /nx; 
dy  =  (yb-ya) /ny; 
dz  =  (zb-za)/nz; 

dt2  =  dt/2; 
dt4  =  dt/4; 

%nxpl  =  nx  +  1; 

%nypl  =  ny  +  1; 

%nzpl  =  nz  +  1; 


denom_obs  =  2  *  delta. ^2; 
denom  =  (2*dt)  *  sigmaa. ^2; 

%const_norm  =  sigmaa...  *  sqrt (2*pi*dt ) ^dim; 

diffus  =  dt2  *  sigmaa.  "^2  /  2; 

diffux  =  diffus (1)  /dx/dx; 

diffuy  =  diffus (2)  /dy/dy; 

diffuz  =  diffus (3)  /dz/dz; 


XX  =  xa  +  dx  *  (0:nx); 
yy  =  ya  +  dy  *  (0 :ny) ; 
zz  =  za  +  dx  *  {0:nz); 


generate.m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%% 

Generate  a  3D 

trajectory  for 

%% 

dX  =  b(X)dt  + 

sigmaa*dW 

%% 

X(0)  =  XO 

%% 

And  generate 

an  observation  for 

%% 

z (k)  =  h(X(t_ 

k) )  +  delta*v(k) 

%% 

Chuanxia  Rao 

%% 

January  1996 

%% 

January  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%  Dependencies: 

%%  initialize .m; 

%%  func_b.m  (function) 

%%  func_h.m  (function) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


% 


%%  Output : 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  initialize; 

disp( '  Generating  the  trajectory  ...'); 

%out_fileZ  =  'data_Z.m'; 

%out_fileX  =  'data_X.m'; 

%Kfinal  =  Kt; 
srdt  =  sqrt(dt); 

X(l, : )  =  XO; 

Z(l, : )  =  ZO; 

randn ( ' seed ' ,  0 )  ; 

V  =  randn (Kf inal , 2 ) ; 

Xk  =  X(l, : ) ; 
for  k  =  l:Kfinal, 

kl=k+l; 

[bkl , bk2 , bk3 ]  =  func_b(Xk(l) ,Xk(2) ,Xk(3) ) ; 
bk  =  [bkl,bk2,bk3] ; 

temp  =  Xk  +  dt*bk;  %%  forward  Euler 

[bkl , bk2 , bk3 ]  =  f unc_b ( temp ( 1 ) ,  temp (2),  temp (3)) 
b2  =  bk  +  [bkl,bk2,bk3] ; 

temp  =  Xk  +  dt2  *  b2 ;  %%  trapezoid 

tempi  =  srdt*diag(sigmaa)*randn(dim,l); 

X(kl,:)  =  temp  +  tempi'; 

Xk  =  X(kl, : ) ; 

[hl,h2]  =  func_h(Xk(l) ,Xk(2) ,Xk(3) ) ; 
temp2  =  diag(delta)*v(k,:)'; 

Z(kl,:)  =  [hl,h2]  +  temp2 ' ; 

end 

%% 

%%  Output  --  graph: 


%% 


C.Rao 

generate.!!! 


kl  =  1 : (Kf inal+1) ; 
kll  =  ( 1 : Kf inal ) +1 ; 
tl  =  0:dt:Tfinal; 
til  =  dt:dt:Tfinal; 

al  =  min(X(kl, 1) ) ;  bl  =  max (X (kl , 1) ) ; 
a2  =  min(X(kl,2) ) ;  b2  =  max (X (kl , 2 ) ) ; 
a3  =  min(X(kl,3) ) ;  b3  =  max (X (kl , 3 ) ) ; 
aZl  =  min(Z (kll, 1) ) ;  bZl  =  max(Z (kll, 1) ) ; 
aZ2  =  min(Z (kll, 2) ) ;  bZ2  =  max (Z (kll , 2 ) ) ; 

f  igure; 

subplot (221) ,  plot(tl,X(kl,l)); 
axis ( [0,Tf inal,  al,bl]); 
subplot (222) ,  plot(tl,X(kl,2)); 
axis ( [0, Tf inal,  a2,b2]); 
subplot (223) ,  plot(tl,X(kl,3)); 
axis ( [0,Tf inal,  a3,b3]); 

subplot (224) ,  plot3 (X(kl, 1) ,  X(kl,2),  X(kl,3),  'r'); 

axis([al,bl,  a2,b2,  a3 , b3 ] ) ; 

f  igure; 

subplot(211)  ,  plot (til, Z (kll, 1) ) ; 
axis ([ 0 , Tf inal,  aZl,bZl]); 
subplot (212) ,  plot (til, Z (kll, 2) ) ; 
axis ( [0, Tf inal,  aZ2,bZ2]); 

clear  temp  tempi  temp2 
clear  bk  bkl  bk2  bk3 
clear  hi  h2 
clear  v  srdt 
clear  al  a2  a3  aZl  aZ2 
clear  bl  b2  b3  bZl  bZ2 
clear  kl  kll 
clear  tl  til 


C.Rao 

moreinit.m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Further  initialization  for  the  new  ADI  for 

%%  dX  =  b(X)dt  +  sigmaa (X) dW,  X(0)  -  N(mO,varO) 

Z(k)  =  h(X(t_k))  +  delta(k)v(k) 

%% - 


%%  Chuanxia  Rao 

%%  May  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Dependencies : 

%%  initialize  .in 

%% - 

%%  Used  in: 

%%  online. m  &  run.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


global  denom_obs 
%global  denom 
%global  nxpl  nypl  nzpl 
global  XX  yy  zz 

global  coefOax  coefObx  coefOcx 
global  coefOay  coefOby  coefOcy 
global  coefOaz  coefObz  coefOcz 
global  coeflax  coeflbx  coeflcx 
global  coeflay  coeflby  coeflcy 
global  coeflaz  coeflbz  coeflcz 

%initialize; 

disp('  Computation  of  the  coefficients  used  in  each  step') 
disp('  and  setting  up  of  the  initial  density  ...') 

[hhl,  hh2]  =  func_h(xx,  yy,  zz) ; 
pk  =  func_pO{xx,  yy,  zz) ; 

%  pk  =  cutsmall (pk) ; 

%  hhl  =  cutsmall (hhl) ; 

%  hh2  =  cutsmall (hh2 ) ; 

%%  hh  =  [hhl,  hh2];  %%  ??? 

%% 

%%  Convection  &  diffusion  coefficients: 

%%  OR  one-dimensional  coefficients: 

%% 


xxl  =  xx(2:nx);  %%  xa  +  dx* ( 1 : nx-1 ) ; 

yyl  =  yy(2:ny);  %%  ya  +  dy*(l:ny-l); 

zzl  =  zz(2:nz);  %%  za  +  dx*(l:nz-l); 

if  (method==01  |  method==02  |  method==ll  |  method==12), 

%%  C~D  splitting  stuff  goes  here  . . . 

disp{'  Warning:  program  not  finished  for  C-D  splitting  ?!') 

elseif  (method==21  |  method==22  |  method==31  |  method==32), 

[  temp_bx ,  temp_by ,  temp_bz  ]  =  f  unc_b  ( xxl ,  yyl ,  z  z  1 )  ; 

[ temp_dbx,  temp_dby ,  temp_dbz ]  =  func_db  (xxl ,  yyl ,  zzl )  ; 

temp_dbx  =  temp_dbx  *  dt2 ; 

temp_dby  =  temp_dby  *  dt2 ; 

temp_dbz  =  temp_dbz  *  dt2 ; 

if  (method==22) ,  %%  central2 

temp_bx  =  temp_bx  *  dt4/dx; 


C.Rao 

moreinit.m 


teinp_by  =  temp_by  *  dt4/dy; 

temp_bz  =  temp_bz  *  dt4/dz; 

teinp_bxp  =  zeros  (size  (temp_bx)  )  ; 

temp_bxn  =  temp_bxp; 

teinp_byp  =  zeros  (size  ( temp_by)  )  ; 

teinp_byn  =  temp^byp; 

temp_bzp  =  zeros (size (temp_bz) ) ; 

temp_bzn  =  temp__bzp; 

elseif  (method==32 ) ,  %%  new2nd 

temp_bx  =  temp_bx  *  dt2/dx; 
temp_by  =  temp_by  dt2/dy; 
temp^bz  =  temp_bz  *  dt2/dz; 
temp_bxp  =  temp_bx  .*  ( temp_bx>0 ) ; 
temp_bxn  =  teinp_bx  .*  ( temp_bx<0 )  ; 
teinp_byp  =  temp_by  .*  ( teinp_by>0 )  ; 
temp__byn  =  temp_by  .*  (teinp_by<0); 
temp_bzp  =  temp^bz  .*  ( temp_bz>0 ) ; 
teinp_bzn  =  temp_bz  (temp_bz<0); 
temp_bx  =  abs ( temp^bx) ; 
temp_by  =  abs ( temp_by) ; 
temp_bz  =  abs  ( teinp_bz )  ; 

elseif  (method==31 ) ,  %%  upwind 

temp_bx  =  temp_bx  *  dt/dx; 
teinp_by  =  temp_by  *  dt/dy; 
temp_bz  =  temp_bz  *  dt/dz; 
temp_bxp  =  temp_bx  .*  ( temp_bx>0 ) ; 
temp_bxn  =  temp_bx  (temp_bx<0) ; 
temp_byp  =  temp_by  .*  (temp_by>0); 
temp_byn  =  temp^by  .*  ( temp_by<0 ) ; 
temp_bzp  =  temp_bz  .*  ( temp_bz>0 ) ; 
teinp_bzn  =  temp_bz  .*  (temp_bz<0)  ; 
teinp_bx  =  abs  ( temp__bx)  ; 
teinp_by  =  abs  { temp__by)  ; 
temp_bz  =  abs  ( teinp_bz )  ; 
temp_dbx  =  temp_dbx  *  2; 
temp_dby  =  temp_dby  *  2 ; 
temp_dbz  =  temp_dbz  *  2 ; 
diffux  =  diffux  *  2; 
diffuy  =  diffuy  *  2; 
diffuz  =  diffuz  *  2; 

elseif  (method==21) ,  %%  centrall 

disp('  Warning:  program  not  finished  for  centrall  ?!'); 

end 


[coef Oax, coef Obx, coef Ocx,  coef lax, coef Ibx, coef lex]  =... 

discret (method,  diffux,  temp_bx, temp_dbx,  temp_bxp,temp_bxn); 
[coef Oay, coef Oby, coef Ocy,  coef lay, coef Iby, coef Icy] 

discret(method,  diffuy,  t emp_by , t emp_dby ,  t emp_byp , t emp_byn ) ; 
[coef Oaz , coef Obz , coef Ocz ,  coef laz , coef Ibz , coef Icz]  =... 

discret (method,  diffuz,  temp__bz , temp_dbz,  temp_bzp, temp_bzn) ; 

else  %%  others 

disp ( '  Discretization  method  not  known  ?!') 

end 

clear  diffus  diffux  diffuy  diffuz 
clear  temp_bx  temp_by  temp_bz 
clear  temp_bxp  temp_byp  temp_bzp 
clear  temp_bxn  temp_byn  temp_bzn 
clear  temp_dbx  temp_dby  temp_dbz 
clear  xxl  yyl  zzl 
clear  dt4  dx  dy  dz 


C.Rao 

discret.m 


function 
%usage : 


[aO,bO,cO,  al,bl,cl] 
[aO,bO,cO,  al,bl,cl] 


discret (method,  dif,  b,  db,  bp,  bn) 
discret (method,  dif,  b,  db,  bp,  bn) 


if (method==22 ) ,  %%  central2 


temp 

=  2*dif  + 

aO  = 

dif 

+  b; 

bO  = 

1  - 

temp; 

cO  = 

dif 

-  b; 

al  = 

-aO ; 

bl  = 

1  + 

temp; 

cl  = 

-cO ; 

elseif (method==32 ) ,  %%  new2nd 

tempi  =  1  +  b; 
temp2  =  2*dif  +  db; 
difn  =  -  dif; 
aO  =  dif  +  bn; 
bO  =  tempi  -  temp2 ; 
cO  =  dif  -  bp; 
al  =  difn  -  bp; 
bl  =  tempi  +  temp2 ; 
cl  =  difn  +  bn; 

elseif (method==31) ,  %%  upwind 

tempi  =  1  +  b; 
temp2  =  2*dif  +  db; 
difn  =  -  dif; 
al  =  difn  -  bp; 
bl  =  tempi  +  temp2 ; 
cl  =  difn  +  bn; 
aO  =  al; 
bO  =  bl; 
cO  =  cl; 


elseif (method==21) , 

disp('  Warning: 


%%  centrall 

program  not  finished  for  centrall  ? ! ' 


else 


%%  others 

disp('  Discretization  method  not  known  ?!') 


end 

return 


C.Rao 

online.m 


a 


•6'0'6-6'0'6'0'6'0'0'0'6'6'0'6'6'0'0'6'0'0'0'6'0'0'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6^'6^'6^^^'6‘6'6'6 


%%  New  ADI  method  for  solving 

%%  dX  =  b(X)dt  +  sigmaa(X)dW,  X{0)  ~  N(mO,varO) 

%%  Z{k)  =  h(X(t_k))  +  delta(k)v(k) 

%% - 

%%  Chuanxia  Rao 

%%  May  1997 

%%  Mar  1997 

%%  Feb  1996 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  %  Dependenc i e  s : 

%%  initialize. m 

%%  moreinit.m 

%%  obs_cor.m 

%%  onestep.m 

%%  cutsmall.m 

%%  generate.m  (for  observations  Z) 

%% - 


%%  Output: 

%%  posterior  density  pkl  (denoted  as  pk) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%' 


%%initialize; 

%%moreinit ; 

%%  --  They  must  be  called  before  running  this  file. 


% 

% 


zl  =  Z (kl, 1) ; 
z2  =  Z (kl, 2) ; 

alpha_kl  =  obs_cor ( zl , z2 ,  hhl,hh2) ; 
alpha_kl  =  cutsmall (alpha_kl) ; 

%%  --  corrector  alpha_kl 

Tpk  =  onestep(pk); 

Tpk  =  cutsmall (Tpk) ; 

%%  --  prior  density  Tpk 

pk  =  alpha_kl  .*  Tpk; 

pk_max  =  max (  max (  max (pk,  [  ] , 3 ) ,  [  ] , 2  ),[],!  ) ; 
i f  pk_max  >  0 , 
pk  =  pk/pk_max; 
end 

pk  =  cutsmall (pk); 


clear 


alpha_kl  Tpk 


%%  --  density  at  step  kl 


C.Rao 
onestep.m 

function  result  =  onestep (method,  u) 

%usage:  result  =  onestep (method,  u) 

%%  Chuanxia  Rao 

%%  Jun  1997 

%%  Mar  1997 

%%  Feb  1996 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  Dependencies: 

%%  moreinit.m 

%%  sol_expl.m 

%%  sol_impl.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if (method  ==  01) , 

%  result  =  split_cd( 

disp ( 'Warning :  Program  for  splitting  not  finished  ?!'); 
elseif (method  ==  02), 

%  result  =  split_dc( 

disp ( 'Warning:  Program  for  splitting  not  finished  ?!'); 
elseif (method  ==  11), 

%  result  =  split_cdc( 

disp ( 'Warning :  Program  for  splitting  not  finished  ? ! ' ) ; 


elseif (method  ==  12), 

%  result  =  split_dcd( 


disp ( 'Warning:  ] 

Program  for  splitting  not 

finished  ? ! ' ) ; 

else 

if (method  ==  22  I  method  ==  32),  %%  second  order  in  time 

global 

coefOax  coefObx 

coef Ocx 

global 

coefOay  coefOby 

coef Ocy 

global 

coefOaz  coefObz 

coef Ocz 

%%  1.1 

(I+AOx)  u: 

u  =  sol_expl(l. 

coef Oax, coef Obx, coef Ocx, 

u)  ; 

%%  1.2 

(I+AOy)  u: 

u  =  sol_expl(2, 

coef Oay, coef Oby, coef Ocy , 

u)  ; 

%%  1.3 

(I+AOz)  u: 

u  =  sol_expl(3, 

coef Oaz , coef Obz , coef Ocz , 

u)  ; 

end 

global 

coeflax  coeflbx 

coef lex 

global 

coeflay  coeflby 

coef Icy 

global 

coeflaz  coeflbz 

coef Icz 

%%  2.3 

(I-Alz) '"{-1}  u: 

u  =  sol_impl(3, 

coeflaz , coeflbz , coef Icz , 

u)  ; 

%%  2.2 

(I-Aly)  "^{-1}  u: 

u  =  sol_impl(2. 

coeflay, coeflby, coef Icy, 

u)  ; 

%%  2.1 

(I-Alx)^{-1}  u: 

u  =  sol_impl(l, 

coeflax, coeflbx, coef lex. 

u)  ; 

result  =  u; 

end 

return 


function 
%  usage: 
%  where 


% 

% 

% 


C.Rao 

sol_expLm 

result  =  sol_expl { i_dir ,  a,  b,  c,  u) 
result  =  sol_expl ( i_dir ,  a,  b,  c,  u) 
i_dir  is  the  direction  to  go  ( l<=i_dir<=3 ) ; 
a,  b,  c  are  three-dimensional  arrays: 
b  is  the  diagonal  and 
u  is  the  right-hand  side. 


%% 

Chuanxia  Rao 

%% 

May  1997 

%% 

Feb  1996 

global 

nx  ny  nz 

%global 

XX  yy  zz 

%global 

time 

%global 

nxpl  nypl  nzpl 

[nxml,  nyml,  nzml]  =  size(b); 
if  (nxml==nx-l  &  nyml==ny-l  &  nzml==nz-l), 
result  =  u; 

%  utemp  =  zeros (b); 

utemp  =  zeros {nxml, nyml, nzml) ; 
if  (i_dir==l), 

for  i=l:nxml, 

utemp(i,:,:)  =a(i,:,:)  .*u(i,  2 : ny , 2 : nz ) . . . 

+  b(i,:,:)  .*  u ( i+1 , 2 : ny, 2 : nz ) . . . 

+  c(i,:,:)  u ( i+2 , 2 :ny, 2 :nz) ; 

end 

elseif  (i_dir==2) , 


for  j=l:nyml 

utemp ( : , j ,  :  ) 

=  a  (  : 

1  j  /  : 

. *  u (2 :nx, j ,  2 :nz) 

+  b(: 

.*  u(2:nx,j+l,2:nz) 

end 

+  c  (  : 

,  j  /  :  1 

1  u (2 :nx, j +2 , 2 :nz ) 

elseif  (i_dir==3), 

for  k=l:nzml, 

utemp (:,:,k)  =  a(:,:,k)  u ( 2 : nx, 2 : ny , k  )... 

+  b(:,:,k)  .*  u ( 2 : nx, 2 : ny , k+1 ) . . . 

+  c(:,:,k)  .*  u (2 :nx, 2 :ny, k+2 ) ; 

end 


else 

end 


disp  (  '  i_dir  is  beyond  1,2,3.') 
result (2 :nx, 2 :ny, 2 :nz)  =  utemp; 


%if  (i_dir==l),  %%  for  nonhomogeneous  boubdary  conditions 

%elseif  (i_dir==2) , 

%elseif  (i_dir==3) , 

%  result (2 :nx, 2 :ny,  1)  =  func_u(time,  xx(2 :nx) ,yy (2 :ny) , zz (1) ) ; 

%  result ( 2 : nx, 2 : ny , nzpl ) =  func_u(time,  xx(2 :nx) ,yy (2 :ny) , zz (nzpl) ) ; 

%end 


else 


end 


disp('  Dimension  problem  in  sol_expl.m') 


function 
%  usage: 
%  where 


C.Rao 

sol_impl.m 

result  =  sol_impl (i_dir,  a,  b,  c,  u) 
result  =  sol_impl (i_dir ,  a,  b,  c,  u) 
i_dir  is  the  direction  to  go  (l<=i_dir<=3 ) ; 
a,  b,  c  are  three-dimensional  arrays: 
b  is  the  diagonal  and 
u  is  the  right-hand  side. 


%%  Chuanxia  Rao 

%%  May  1997 

%%  Feb  1996 


global  nx  ny  nz 
%global  nxpl  nypl  nzpl 

result  =  u; 

utemp  =  u ( 2 : nx , 2 : ny , 2 : nz ) ; 

%if  (i_dir==l) ,  %%  for  nonhomogeneous  boubdary  conditions 

%elseif  (i_dir==2), 

%elseif  (i_dir==3), 

%  utemp =  utemp -  a(:,:,  1)  .*  u ( 2 : nx, 2 : ny , 1 ) ; 

%  utemp nzml )  =  utemp nzml) -  c(:,:,nzml)  .*  u (2 :nx, 2 :ny, nzpl) 

/ 

%end 

utemp  =  tridiag3d (i_dir ,  a,b,c,  utemp); 
result (2 :nx, 2 :ny, 2 :nz)  =  utemp; 


C.Rao 
tridiagSd.m 

result  =  tridiag3d(i_dir,  a,  b,  c,  d) 
result  =  tridiag3d (i_dir ,  a,  b,  c,  d) 
i_dir  is  the  direction  to  go  (l<=i_dir<=3 ) ; 
a,  b,  c,  d  are  three-dimensional  arrays  of  the  same  size 
b  is  the  diagonal,  and  d  is  the  right-hand  side. 

(d  is  changed  after  the  call)  ?  -  No. 

%%  Chuanxia  Rao 

%%  May  1997 

%%  Feb  1996 

%global  dim 
%dim  =  3 ; 

na  =  size (a) ; 
nb  =  size (b) ; 
nc  =  size (c) ; 
nd  =  size (d) ; 

if(na==nb  &  nb==nc  &  nc==nd  &  length (nd) <=3 ) 
n  =  na(i_dir); 
if (i_dir==l) 

V  ( 1 ,  ;  ,  :  )  =  b  ( 1 ,  :  ,  :  )  ; 

for  i=2:n, 

xmult  =  a(i,:,:)  ./  v(i-l,:,:); 

v(i,:,;)  =b(i,:,:)  -  xmult  .*  c(i-l,:,;); 
d(i,:,:)  =d(i,:,:)  -  xmult  .*  d(i-l,:,:); 
end 

result (n, :,; )  =  d(n,:,:)  ./  v(n, :,:); 

for  i=(n-l) 

d(i,:,:)  =d(i,:,:)  -c(i,:,:)  .*  result(i+l,:,:); 

result (i, :,: )  =d(i,:,:)  ./v(i,:,:); 

end 

elseif ( i_dir==2 ) 

v( : ,1, : )  =  b( : ,1, : ) ; 
for  i=2:n, 

xmult  =  a(:,i,:)  ./  v(:,i-l,:); 

v(:,i,:)  =b(:,i,:)  -  xmult  .*  c(:,i-l,:); 
d(:,i,:)  =d(:,i,:)  -  xmult  .*  d(:,i-l,:); 
end 

result (:, n, : )  =  d(:,n,:)  ./  v(:,n,:); 

for  i= (n-1) ; -1 : 1 , 

d(:,i,:)  =  d(:,i,:)  -  c(:,i,:)  .*  result(:,i+l,:); 

result (;, i, : )  =d(:,i,:)  ./v(:,i,:); 

end 

elseif { i_dir==3 ) 

v(:, :,1)  =  b(:, :,1); 
for  i=2:n, 

xmult  =  a(:,;,i)  ./  v(:,;,i-l); 

v(:,:,i)  =b(:,:,i)  -  xmult  .*  c(:,:,i-l); 
d(:,:,i)  =d(:,:,i)  -  xmult  .*  d(:,:,i-l); 
end 

result  ,n)  =  d(:,;,n)  ./  v(:,:,n); 


function 
%  usage : 
%  where 
% 

% 

% 


C.Rao 

tridiagSd.m 


else 

end 

elseif ( 

else 

end 


for  i= (n-1) : -1 : 1, 

d(:,:,i)  =d(:,:,i)  .*  result(:,:,i+l); 

result i)  =d(:,:,i) 

end 

disp ( '  i_dir  is  beyond  1,2,3. ' ) 

length (nd) >3  ) 
disp('  dimension  >  3') 

disp('  Dimensions  do  not  agree  in  tridiag3d.') 


C.Rao 
obs_cor.m 

function  result  =  obs_cor (zl, z2,  hhl,hh2) 

%  usage:  result  =  obs_cor ( zl , z2 ,  hhl,hh2) 

%  usage:  result  =  obs_cor(z,  hh) 

%  where  z  is  the  measurement,  and 
%  hh  is  the  part  without  noise, 

dim  =  3 ;  dim_obs  =  2 ; 

Chuanxia  Rao 
May  1997 
Feb  1996 

Used  in:  online. m 
global  denom_obs  %%  defined  in  moreinit.m 

temp  =  zl  -  hhl 

tempi  =  temp/denom_obs (1)  .*  temp; 

temp  =  z2  -  hh2 

tempi  =  tempi  +  temp/denom_obs (2 )  .*  temp; 

tempi  =  tempi  -  min {  min {  min ( tempi ,[],3),[],2  ),[],! 

result  =  exp (  -  tempi  ); 


function 
%usage : 

% 

% 


C.Rao 

plotting.m 

plotting (xyz,  pk) 
plotting (xyz,  pk) 
where  length (xyz)  =  dim (=3) 

size(pk)  =  [nxpl  nypl  nzpl] 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%  This  is  for  plotting  the  tracking  process. 

%%  Chuanxia  Rao 

%%  April  1997 

%%  May  1997 

%% 


global  xa  ya  za 
global  xb  yb  zb 
global  XX  yy  zz 
global  resol_vert 

%figure 


X  =  xyz ( 1 ) ; 
y  =  xyz ( 2 ) ; 
z  =  xyz ( 3 )  ; 
pk_xy  =  sum(pk,  3); 
pk_x  =  sum(pk_xy,  2); 

subplot (221) 

plot(xx,  pk_x) ; 
hold  on; 

lower  =  0 ; 

upper  =  max(pk_x) ; 

upper  =  upper  *  1.1; 

resol  =  resol_vert  *  upper; 

vert  =  lower:resol:upper; 

plot (x*ones (size (vert) ) ,  vert,  'g'); 

axis ( [xa  xb  lower  upper] ) ; 

pause (1) ; 

hold  off; 

subplot (222 ) 

pk_xy  =  pk_xy' ; 
mesh(xx,  yy,  pk_xy) ; 
hold  on; 

m  =  max (max (pk_xy) ) ; 
m  =  m  *  1.1; 

plot3 ( [x  x] ,  [y  y] ,  [0  m] ,  'g' ,  'LineWidth' , [2] ) ; 

axis ( [xa  xb  ya  yb  0  m] ) ; 
pause ( 1 ) ; 
hold  off; 


C.Rao 
plotting.m 

subplot (224) 

contour (xx,  yy,  pk_xy) ; 
hold  on; 

plot(x,  y,  'g*' ) ; 
pause ( 1 ) ; 
hold  off; 

subplot  (223 ) 

plots (x,  y,  z,  'g*'); 
axis ( [xa  xb  ya  yb  za  zb] ) ; 
hold  on; 

[one,  ijk]  =  max (pk, [ ] , 3 )  ; 

[ one ,  i j ]  =  max ( one , [ ] , 2 ) ; 

[ one ,  i ]  =  max ( one ,  [  ]  ,  1 )  ; 

X  =  XX (i) ; 

y  =  yy(ij (i) ) ; 

z  =  zz (ijk(i, ij (i) )  )  ; 

%  plots (x,  y,  z,  'ro',  ' LineWidth' , [2 ] ) ; 

plots (x,  y,  z,  ' r* ' ) ; 
grid  on; 
pause ( 1 ) ; 

%  hold  off; 


return 


C.Rao 

cutsmall.m 

function  result 

=  cutsmall (array) 

%  usage 

:  result 

=  cutsmall (array) 

%% 

Used  in: 

online .m 

%% 

Chuanxia 

Rao 

9'9' 

May  1997 

global 

small 

%%  defined  in  initial 

result  = 

(array  >  small)  .*  array; 

% 

result  = 

sparse (  result  ) ; 

ize  ,m 


C.Rao 

advection.m 

function  [x_down, Y_down, z_down]  =  advection{dt. 

%  usage 

:  [ x_down , y_down , z_down ]  =  advec  t i on ( dt , 

%% 

Chuanxia  Rao 

%% 

Februry  1997 

%% 

April  1997 

%% 

June  1997 

global 

substeps 

-  Y_up ,  z_up ) 
X— up ,  y_up ,  z_up ) 


%dt6  =  dt/6; 
dt2  =  dt/2; 

X  =  x_up; 
y  =  x_up; 
z  =  x_up ; 


for  k 

=  1: substeps. 

[bkl,bk2,bk3]  =  func_b (x,y, z) 

/ 

[bl,b2,b3]  =  func_b (x+dt*bkl, 
X  =  X  +  dt2  *  (bkl  +  bl) ; 
y  =  y  +  dt2  *  (bk2  +  b2 )  ; 

y+dt* 

bk2 , z+dt*bk3 ) ; 

z  =  z  +  dt2  *  (bk3  +  b3 ) ; 

%% 

trapezoid 

% 

bk  =  func_b(x); 

% 

temp  =  func_b(x  +  dt2*bk); 

% 

X  =  X  +  dt  *  temp; 

%% 

modified  Euler 

%% 

kl  =  func_b(x); 

%% 

k2  =  func_b(x  +  dt2*kl) ; 

%% 

k3  =  func_b(x  +  dt2*k2); 

%% 

k4  =  func_b(x  +  dt*k3); 

%% 

temp  =  kl  +  2*k2  +  2*k3  +  k4; 

%% 

X  =  X  +  dt6  *  temp; 

%% 

Runge-Kutta-4 

end 

x_down  =  X; 
y_down  =  y; 
z_down  =  z; 

C.Rao 
func_b.m 

function  [bl,b2,b3]  =  func_b(x,y, z) 

%  usage:  [bl,b2,b3]  =  func_b (x, y , z ) 

%  where  x,y,z  are  vectors. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  signal  propogation  function 

%%  Chuanxia  Rao 

%%  Februry  1997 

%%  May  1997 

%% - 

%%  Called  in: 

%%  moreinit.m 

%%  generate,  m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


for  j=l : length (y) , 
for  k=l : length (z) , 

%  bl  (1:  length  (x)  ,  j,  k)  =  200*y(j)''3  + 

bl  (1 :  length  (x)  ,  j,  k)  =  -  200*y(j)''3 

end 

end 

bl  =  bl/2; 


% 

% 


b2  =  1./2  *  ones (size (bl) ) ; 
b2  =  -1./2  *  ones (size (bl) ) ; 

b3  =  1./4  *  ones (size (bl) ) ; 
b3  =  -1./4  *  ones (size (bl)); 


b  = 


50*z(k)  -  50; 
+  50*z (k) ; 


[bl,  b2,  b3]; 


C.Rao 
func_db.m 

function  [dbl, db2 , db3]  =  func_db(x,y, z) 
%  usage:  [dbl , db2 , db3 ]  =  func_db (x, y, z ) 
%  where  x,y,z  are  vectors. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  derivatives  dxbl,  dyb2 ,  dzb3 

%%  (The  scalar  function  Delta  *  b) 

%%  Chuanxia  Rao 

%%  March  1997 

%%  May  1997 


%%  Called  in: 

%%  moreinit .m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


nl 

n2 

n3 


%n 


% 

% 


length (x) ; 
length (y) ; 
length ( z ) ; 

dbl  =  zeros (nl , n2 , n3 ) ; 
db2  =  zeros (nl,n2,n3) ; 
db3  =  zeros (nl,n2,n3) ; 
length(x)  *  length(y)  *  length(z); 
dbl  =  zeros (n,l); 
db2  =  zeros (n,l); 
db3  =  zeros (n,l); 
db  =  zeros (n, 1) ; 


%% 

db  : 

=  dbl  + 

db2  +  db3; 

%% 

dbl 

=  db  / 

3; 

%% 

db2 

=  dbl; 

%% 

db3 

=  dbl; 

C.Rao 
func_h.m 

function  [hl,h2]  =  func_h(x,y, z) 

%  usage;  [hl,h2]  =  func_h (x, y , z ) 

%  where  x,y,z  are  vectors. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  observation  function 

%%  Chuanxia  Rao 

%%  April  1997 

%%  May  1997 

%% - 


%%  Called  in: 
%%  moreinit.m 
%%  generate,  m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%  X  =  [x(:,l),  x(:,2),  x(:,3)] 

for  i=l : length (x) , 
for  j=l : length (y) , 

templ2  =  x(i)''2  +  y(j)''2; 
if  templ2==0, 
for  k=l ; length ( z ) , 

hl(i,j,k)  =  0; 

h2(i,j,k)  =  asin(  sign(z(k))  ); 

end 

else 

argl  =  x{i)  /  sqrt  ( teinpl2 )  ; 
tempi  =  sign(y(j))  *  acos(argl); 
for  k=l : length ( z ) , 

hl(i,j,k)  =  tempi; 
temp  =  templ2  +  z(k)^2; 
arg2  =  z(k)  /  sqrt(temp); 
h2(i,j,k)  =  asin(arg2); 

end 

end 

end 

end 

%  argl  =  x(:,l)  ./  sqrt  (x  (  :  ,  1)  .  “^2  +  x(:,2).^2); 

%  arg2  =  x(:,3)  ./  sqrt  (x  (  :  ,  1)  .  ^^2  +  x(:,2).''2  +  x(:,3).''2); 

%if(x(:,  2X0), 

%%  hi  =  2*pi  -  acos{argl); 

%  hi  =  acos(argl); 

%else 

%  hi  =  acos(argl); 

%end 

%  h2  =  temp  .  *  h2  ; 

%  h  =  [hi,  h2]  ; 


function 
%  usage : 
%  where 


pO  =  func_pO (x, y, z ) 
pO  =  func_pO (x, y, z ) 
x,y,z  are  vectors. 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%  The  initial  density  function 

%%  Chuanxia  Rao 

%%  Februry  1997 

%%  May  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Called  in: 


%%  moreinit.m 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


global  mO  varO_inv 
%  temp  =  sqrt (2*pi) ^3 ; 

%  tempo  =  temp*sqrt (varO_inv(l) *var0_inv(2) *var0_inv(3 ) ) ; 

for  i=l : length (x) , 

tempi  =  (x(i) -mO  (1)  ) '^2  *  varO_inv(l)  ; 
for  j=l : length (y) , 

temp2  =  (y ( j ) -mO (2) ) ^2  *  var0_inv(2)  ; 
for  k=l : length (z)  , 

temp3  =  (z  (k) -mO  (3  )  ) ''2  *  var0_inv(3)  ; 
temp3  =  tempi  +  temp2  +  temp3 ; 
pO(i,j,k)  =  exp (  -temp3/2  ); 

end 

end 

end 

%  pO  =  pO  /  tempo ; 

pO  =  pO  /  max (  max (  max (pO , [ ] , 3 ) , [ ] , 2  ),[],!  ) ; 


function  b  =  func_b(x) 

%  usage:  b  =  func_b(x) 

%  where  x  is  a  matrix  of  size  (:,  dim=3). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  signal  propogation  function 
%%  Chuanxia  Rao 

%%  Februry  1997 

%%  April  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Called  in: 

%%  advection.m  (function) 

%%  generate.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  x  =  [x(:,l),x(:,2),x(:,3)] 

global  gam 

bl=-x(:,2); 

b2  =  -exp(-gam*x(:,l))  .*  x(:,2).^2  .*  x(:,3); 
b3  =  zeros(size(bl)); 


b  =  [bl,b2,  b3]; 
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function  db  =  func_db(x) 

%  usage:  db  =  func_db(x) 

%  where  x  is  a  matrix  of  size  (:,  dim=3). 

%  db  is  a  vector  of  size  (:). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  scalar  function  Delta  *  b 
%%  Chuanxia  Rao 

%%  April  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Called  in: 

%%  split_cdc.m 

%%  run.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  x  =  [x(:,l),x(:,2),x(:,3)] 

global  gam 

%  dbl  =  zeros(size(x(:,l)); 

%  db2  =  -2*  exp(-gam*x(:,l))  .*  x(:,2)  .*  x(:,3); 

%  db3  =  zeros(size(x(:,3)); 

%  db  =  dbl +db2  +  db3; 

db  =  -2  *  exp(-gam*x(:,l))  .*  x(:,2)  .*  x(:,3); 
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function  h  =  func_h(x) 

%  usage:  h  =  func_h(x) 

%  where  x  is  a  matrix  of  size  (:,  dim=3). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  The  observation  function 

%%  Chuanxia  Rao 

%%  April  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%% 

Called  in: 

%% 

offline.m 

%% 

generate.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  x  =  [x(:,l),x(:,2),x(:,3)] 
global  MMH 

temp  =  MM  +  (x(:,l)  -  H)A2; 
h  =  sqrt(temp); 
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function  pO  =  func_pO(x) 

%  usage:  pO  =  func_pO(x) 

%  where  x  is  x  is  a  matrix  of  size  (:,  dim=3). 

%  pO  is  a  vector  of  size  (:). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%  The  initial  density  function 
%%  Chuanxia  Rao 

%%  Februry  1997 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%  Called  in: 

%%  offline.m 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
global  mO  varO_inv 

temp  =  sqrt(2*pi)''3; 

temp  =  temp*sqrt(varO_inv(l)*varO_inv(2)*varO_inv(3)); 

tempi  =  (x(:,l)-m0(l)).'^2  *  varO_inv(l); 

temp2  =  (x(:,2)-m0(2)).''2  *  var0_inv(2); 

temp3  =  (x(:,3)-m0(3)).'^2  *  var0_inv(3); 

temp3  =  tempi  +  temp2  +  temp3; 

pO  =  exp(  -temp3/2  )  /  temp; 

pO  =  pO  /  max(pO); 
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! !  Modules  for  initial  data: 

I  I 

MODULE  init_datal 
REAL*  8  xi 
INTEGER  ix 

INTEGER,  PARAMETER  : :  dimen  =  1 
INTEGER,  PARAMETER  : :  nx  =  64 
INTEGER,  PARAMETER  : :  nxl  =  nx-1 
!  INTEGER,  PARAMETER  : :  nx2  =  nx-2 

INTEGER,  PARAMETER  : :  mt  =  256 
END  MODULE  init_datal 

MODULE  init_data2 
USE  init_datal 
REAL* 8  dt,  dx 
REAL *8  dt2,  dtx2 
REAL* 8  to,  tl 

REAL* 8  xa,  xb,  ya,  yb,  za,  zb 
REAL*8  xx{0:nx) 

PARAMETER  (tO=ODO,  tl=lD0) 

PARAMETER  (xa=-1.0DO,  xb=1.0D0) 

PARAMETER  (dx  =  (xb-xa) /nx) 

!  PARAMETER  (dt  =  IDO/mt) 

PARAMETER  (dt  =  (tl-tO)/mt) 

PARAMETER  (dt2  =  dt/2) 

PARAMETER  (dtx2  =  dt2/dx) 

PARAMETER  (xx=xa+dx*  (/(i,  i=0,nx)/)  ) 

END  MODULE  init_data2 

MODULE  init_data3 
USE  init_data2 

REAL* 8  diffus,  ddiff,  dd_neg 
REAL* 8  dneg,  dpos 
PARAMETER  (  diffus  =  IDO  ) 

PARAMETER  (  ddiff  =  dif fus*dtx2 /dx  ) 
PARAMETER  (  dd_neg  =  -  ddiff  ) 

PARAMETER  (  dneg  =  IDO  -  2  *  ddiff  ) 
PARAMETER  (  dpos  =  IDO  +  2  *  ddiff  ) 

END  MODULE  init_data3 


MODULE  coeffunc 

REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
CONTAINS 

FUNCTION  func_b(x) 

!  1 

!  !  The  coefficient  functions  for  convection  (or  drift)  term 

[  I 

REAL*8  func_b,  x 
!  REAL* 8  X 

!  REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
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INTRINSIC  DSIN 
func_b  =  -3D0  *  DSIN (PI  *  x) 
func_b  =  ODO 
RETURN 

END  FUNCTION  func_b 
FUNCTION  func_c{x) 


The  coefficient  function  for  u  (zero-th  order  term) 

REAL*8,  PARAMETER  : :  PI  =  3,14159265358979323 
REAL *8  X,  temp 
REAL* 8  func_c,  x,  temp 
INTRINSIC  DCOS 
temp  =  3D0  *  DCOS (PI  *  x) 
temp  =  temp  -  PI 
func_c  =  PI  *  temp 
func_c  =  -PI  *  PI 
func_c  =  0 
RETURN 

END  FUNCTION  func_c 

FUNCTION  func_f ( t,x,y, z) 

The  force  function  f: 

REAL*  8  func_f,  t,  x,  y,  z 

RETURN 

END  FUNCTION  func_f 
FUNCTION  f unc_u ( t ,  x ) 


The  true  solution  u(t,  x) : 

REAL *8  func_u,  t,  x 
REAL* 8  t,  X 

REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
INTRINSIC  DSIN,  DEXP 
func_u  =  DSIN{PI*x)  *  DEXP ( -2 *PI*PI* t ) 
func_u  =  DSIN(PI*x)  *  DEXP ( -PI*PI*t ) 

RETURN 

END  FUNCTION  func_u 
END  MODULE  coeffunc 

PROGRAM  adild 

Main  program  begins : 

USE  init_data2 
USE  coeffunc 
INTEGER  Kmax 
REAL *8  time 
REAL* 8  ul(0:nx) 

REAL*8  u2(0:nx) 

REAL *8  t_fin,  v 
REAL* 8  errorl,  error2 ,  value 
REAL* 8  ABS,  MAX 

REAL  cputime,  tarray(2),  ETIME  !,  DTIME 
EXTERNAL  ETIME  ! ,  DTIME 
EXTERNAL  func_u 
EXTERNAL  march 
EXTERNAL  solution 
INTRINSIC  DABS,  DMAXl ,  IDINT 
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func_u{t,x)  =  DSIN(PI*x)  *  DEXP(-PI*PI*t) 
func_u(t,x)  =  DSIN(PI*x)  *  DEXP(-2*PI*PI*t) 

cputime  =  ETIME { tarray) 
cputime  =  DTIME (tarray) 


Kmax  =  2 
t_fin  =  Kmax*dt 
t_fin  =  tl 

CALL  solution ( t_f in,  ul ,  u2 ) 

!!  "t_fin''  might  be  slightly  changed  after  the  CALL. 

DO  ix  =  0 ,  nx 
xi  =  XX (ix) 

ul(ix)  =  func_u(tO,  xi) 
u2(ix)  =  func_u(tO,  xi) 

END  DO 


New  method: 

CALL  march (  Kmax ,  1 ,  ul  ) 


Old  method: 

CALL  march (  Kmax,  2 ,  u2  ) 

cputime  =  DTIME ( tarray) 
cputime  =  ETIME ( tarray)  -  cputime 


Maximum  errors : 

value  =  ODO 
errorl  =  ODO 
error2  =  ODO 
DO  ix  =  1,  nxl 
xi  =  xx(ix) 

V  =  func_u ( t_f in,  xi) 
value  =  DMAXl (  value,  DABS (v)  ) 

errorl  =  DMAXl (  errorl,  DABS (  v  -  ul(ix)  )  ) 

error2  =  DMAXl (  error2 ,  DABS (  v  -  u2(ix)  )  ) 
END  DO 

PRINT  * ,  '  dx  =  ' ,  dx 

PRINT  * ,  '  dt  =  ' ,  dt 
PRINT  '  tl  =  ',  tl 

PRINT  *,  '  t_fin  =  t_fin 

PRINT  *,  '  total  cpu  =  cputime 

PRINT  *,  '  max  value  =  value 

PRINT  * ,  '  max  error_new  =  ' ,  errorl 

PRINT  * ,  '  max  error_old  =  ' ,  error2 

END  PROGRAM  adild 

SUBROUTINE  march (Kmax,  label,  u) 

USE  coefficients 
USE  coeffunc 
USE  init_data3 
INTEGER  label 
INTEGER  Kmax,  ktime 
LOGICAL  judge 
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REAL* 8  temp_b,  terap_c,  vect_b 
REAL  *8  time 
REAL*8  u(0:nx) 

REAL *8  ux(nxl) 

REAL*8  coefOa(nxl) 

REAL*8  coefOb{nxl) 

REAL*8  coefOc(nxl) 

REAL*8  coefla{nxl) 

REAL *8  coeflb(nxl) 

REAL*8  coef Ic (nxl ) 

REAL*8  aax(nxl) ,  bbx{nxl) ,  ccx{nxl) ,  ddx(nxl) 

REAL  *8  coefld(nxl) 

REAL* 8  f{0:nx) 

EXTERNAL  func_u 
EXTERNAL  f unc_b ,  f unc_c 
INTRINSIC  DABS 
EXTERNAL  coef 
INTENT (IN)  ::  Kmax,  label 
INTENT (INOUT)  ::  u 

func_b(x)  =  -3D0  *  DSIN(PI  *  x) 
func_b(x)  =  ODO 

func_c(x)  =  PI  *  (3DO*DCOS (PI*x)  -  PI) 
func_c(x)  =  -PI  *  PI 
func_c(x)  =  ODO 

CALL  coef (label,  coef Oa , coef Ob, coef Oc ,  coef la, coef lb, coef Ic) 
CALL  coef (label) 

DO  ix  =  1,  nxl 
xi  =  XX (ix) 

vect_b  =  func_b(xi)  *  dtx2 
temp_c  =  func_c(xi)  *  dt2 
judge  =  (vect_b  .LE,  ODO) 
t  emp_b  =  DAB  S ( ve c  t _b ) 

IF  (label==l)  THEN 
IF  (judge)  THEN 

coefOa(ix)  =  ddiff 
coefOc(ix)  =  ddiff  +  vect__b 
coefla{ix)  =  dd_neg  +  vect_b 
coeflc(ix)  =  dd_neg 
ELSE 

coefOa(ix)  =  ddiff  -  vect_b 
coefOc(ix)  =  ddiff 
coefla(ix)  =  dd_neg 
coeflc(ix)  =  dd_neg  -  vect_b 
END  IF 

coefOb(ix)  =  dneg  +  temp_b  +  temp_c 

coeflb(ix)  =  dpos  +  temp_b  -  temp_c 

ELSE 

vect_b  =  vect_b  /  2 
coefOa{ix)  =  ddiff  -  vect_b 

coefOc(ix)  =  ddiff  +  vect_b 

coefla(ix)  =  dd_neg  +  vect_b 
coeflc(ix)  =  dd_neg  -  vect_b 
coefOb(ix)  =  dneg  +  temp_c 
coeflb(ix)  =  dpos  -  temp_c 

END  IF 
END  DO 

DO  ktime  =  1,  Kmax 

time  =  to  +  dt*ktime 
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& 

& 


DO  ix  =  1 , nxl 

aax{ix)  =  coefla(ix) 
bbx(ix)  =  coeflb{ix) 
ccx{ix)  =  coeflc(ix) 
ddx(ix)  =  coefOa{ix) 
+  coefOb(ix) 
+  coefOc(ix) 

END  DO 


*  u(ix-l) 

*  u{ix) 

*  u(ix+l) 


u{l:nxl)  =  ux 

u ( 0 )  =  f unc_u ( t ime ,  xx ( 0 ) ) 
u(nx)  =  func_u{time,  xx(nx)) 
u  =  u  +  f (ktime,  : ) 

+  f (ktime,  ix) 


& 

Sc 


aax  =  coef la (2 :nxl ) 
bbx  =  coef lb (1 : nxl) 
ccx  =  coeflc (1 :nx2 ) 
ddx  =  u { 1 : nxl ) 

ddx(l)  =  ddx(l)  -  coeflad,  1,  jy,kz)  *u(0,  jy,kz) 
ddx{nxl)  =  ddx(nxl)  -  coeflcd,  nxl ,  jy,  kz )  *u  (nx,  jy,  kz) 
CALL  tridiag(nxl,  aax,  bbx,  ccx,  ddx,  ux) 
ud:nxl)  =  ux 


END  DO 
RETURN 

END  SUBROUTINE  march 
SUBROUTINE  tridiag(N,  A,B,C,D,  X) 


The  tridiagonal  solver: 

INTEGER  N 

REALMS  A(N)  ,B(N)  ,C(N)  ,D(N)  ,X(N) 

REAL*8  A(N-l)  ,B{N)  ,C(N-1)  ,  D(N),  X(N) 
REAL *8  XMULT 
INTENT (IN)  ::  N,  A,C 
INTENTdNOUT)  ::  B,D,  X 
DO  2  I  =  2,N 

XMULT  =  A(I) /B(I-l) 

XMULT  =  A(I-l) /B(I-l) 

Bd)  =  B(I)  -  XMULT*C(I-1) 

D(I)  =  D(I)  -  XMULT*D(I-1) 

2  CONTINUE 

X(N)  =  D(N) /B(N) 

DO  3  I  =  N~l,l,-1 

X(I)  =  (D{I)  -  C(I)*X(I+1) )/B(I) 

3  CONTINUE 
RETURN 

END  SUBROUTINE  tridiag 


C.Rao 
adi3D.f90 

I  !  !  ]  1  !  !  i  I  !  I  !  !  !  !  j  I  I  1  I  I  I  I  1  1  I  r  I  I  I  I  j  I  t  i  !  1  I  i  j  I  I  t  j  I  I  I  J  t  I  j  j  I  j 
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Modules  for  initial  data: 

MODULE  init_datal 

REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
REAL*8  xi ,  y j ,  zk 
INTEGER  ix,  jy,  kz 
INTEGER,  PARAMETER  : :  dimen  =  3 
INTEGER,  PARAMETER  ::  nx=32,  ny=32,  nz=32 
INTEGER,  PARAMETER  ::  nx=64,  ny=64,  nz=64 
INTEGER,  PARAMETER  ::  nxl=nx-l,  nyl=ny-l,  nzl=nz-l 
INTEGER,  PARAMETER  ::  nx2=nx-2,  ny2=ny-2,  nz2=nz-2 
INTEGER,  PARAMETER  : :  mt  =  100 
INTEGER,  PARAMETER  : :  mt  =  200 
INTEGER,  PARAMETER  : :  msubt  =  4 
END  MODULE  init_datal 

MODULE  init_data2 

USE  init_datal 

REAL  *  8  di f  f us ( dimen ) 

REAL* 8  dt,  dx (dimen) 

REAL*8  dtdim2,  dtx2 (dimen) 

REAL* 8  to,  tl 

REAL *8  xa,  xb,  ya,  yb,  za,  zb 
REAL*8  xx(0:nx) 

REAL* 8  yy(0:ny) 

REAL*8  zz(0:nz) 

PARAMETER  (t0=0D0,  tl=lD0) 

PARAMETER  (xa=-1.0D0,  xb=1.0D0) 

PARAMETER  (ya=-1.0D0,  yb=1.0D0) 

PARAMETER  (za=-1.0D0,  zb=1.0D0) 

PARAMETER  (  dx  =  (/  (xb-xa) /nx,  (yb-ya) /ny,  (zb-za)/nz  /)  ) 

PARAMETER  (dt  =  iDO/mt) 

PARAMETER  (dt  =  (tl-t0)/mt) 

PARAMETER  (dtdim2  =  dt/dimen/2) 

PARAMETER  (  dtx2  =  dt/dx/2  ) 

PARAMETER  (subdt  =  dt /msubt)  ! ! (subdt  =  dt) 

PARAMETER  (xx=xa+dx(l)  *  (/(i,  i=0,nx)/)  ) 

PARAMETER  (yy=ya+dx(2)  *  (/(j,  j=0,ny)/)  ) 

PARAMETER  (  zz  =  za  +  dx(3)  *  (/(k,  k=0,nz)/)  ) 

PARAMETER  (  diffus  =  (/  1D0/6D0,  1D0/6D0,  1D0/6D0  /)  ) 
PARAMETER  (  diffus  =  (/  IDO,  IDO,  IDO  /)  ) 

DATA  diffus  /IDO,  2D0,  3D0/  !!  diffusion  parameters 
END  MODULE  init_data2 

MODULE  init_out 

CHARACTER  (*)  ::  file_out  !,  fmt_out,  status__out 

PARAMETER  (file_out  =  ' result . temp ' ) 

END  MODULE  init_out 


MODULE  coeffunc 
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USE  init_datal 
CONTAINS 

The  coefficient  functions  for  convection  (or  drift)  term: 

FUNCTION  func_b{x,y, z) 

REAL  *  8  f unc_b ( dimen ) 

REAL* 8  X,  y,  z 
func_b(l)  =  -  DSIN(PI  *  x) 
func_b(2)  =  -  DSIN(PI  *  y) 
func_b{3)  =  -  DSIN{PI  *  z) 

END  FUNCTION  func_b 

The  coefficient  function  for  u  (zero-th  order  term) : 

FUNCTION  func_c (x, y , z ) 

REAL* 8  func_c,  x,  y,  z,  temp 

temp  =  DCOS{PI  *  x) 
temp  =  temp  +  DCOS(PI  *  y) 

temp  =  temp  +  DCOS(PI  *  z) 

temp  =  temp  ~  PI 

func_c  =  PI  *  temp 

END  FUNCTION  func_c 

The  force  function  f: 

FUNCTION  func_f {t,x,y,z) 

REAL* 8  func_f,  t,  x,  y,  z,  temp 
END  FUNCTION  func_f 

FUNCTION  func__u{t,  x,y,z) 

REAL *8  func_u,  t,  x,  y,  z,  temp 
REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
INTRINSIC  DSIN,  DEXP 
temp  =  DSINCPI  *  x) 
temp  =  temp  *  DSIN(PI  *  y) 
temp  =  temp  *  DSIN(PI  *  z) 
func_u  =  temp  *  DEXP (-3D0/2*PI*PI* t ) 
func_u  =  temp  *  DEXP ( -4*PI*PI*t ) 

END  FUNCTION  func_u 
END  MODULE  coeffunc 

PROGRAM  adi3d 

Main  program  begins : 

USE  init_out 
USE  init_data2 
USE  coeffunc 

REAL* 8  u { 0 :nx, 0 :ny, 0 :nz) 

REAL* 8  t_fin,  v,  error,  value 
REAL *8  ABS,  MAX 
REAL  cputime,  tarray{2),  ETIME 
EXTERNAL  ETIME 
EXTERNAL  solution 
INTRINSIC  DABS,  DMAXl 

INTRINSIC  ABS,  EXP,  SQRT,  ATAN,  IDINT,  MAX,  MATMUL, 

!  MAXLOC,  DOT^PRODUCT 

PRINT  *,  'Off-line  calculations,  please  wait 

OPEN  (UNIT=8,  FILE=file_out,  STATUS= ' REPLACE ' ) 

WRITE  (8,  *)  'nx  =  ' ,  nx,  ' ; ' 

WRITE  (8,  *)  'ny  =  ',  ny, 
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WRITE  (8,  *)  'nz  =  nz, 

WRITE  (8,  *)  'mt  =  mt, 

!  WRITE  (8,  '(9A)')  ' index  =  [' 

!  333  FORMAT  (A4,  14,  Al) 

cputime  =  ETIME ( t array) 

!  cputime  =  DTIME ( tarray) 


t_fin  =  tl 
t_fin  =  2*dt 
CALL  solution {t_f in,  u) 

! !  "t_fin"  might  be  slightly  changed  after  the  CALL. 


cputime  =  DTIME ( tarray) 
cputime  =  ETIME ( tarray)  -  cputime 


PRINT 

* 

t 

'  dx  =  ' 

' ,  dx 

PRINT 

★ 

'  dt  =  ^ 

'  ,  dt 

PRINT 

★ 

/ 

'  tl  =  ' 

tl 

PRINT 

* 

/ 

'  t_fin 

..fin 

PRINT 

* 

'  total 

cpu  =  ' 

'  ,  cputime 

WRITE  (8,  11)  'dx=',  dx(l), 


WRITE 

(8, 

*) 

'dx  = 

dx{l) 

f  »  t 
!  / 

WRITE 

(8, 

*) 

'  dy  =  ' 

dx(2) 

t  .  f 
/  / 

WRITE 

(8, 

*) 

'dz  =  ' 

dx(3) 

f  »  f 
/  / 

WRITE 

(8, 

*) 

'dt  =  - 

dt,  ' 

,  / 
t 

WRITE 

(8, 

,  22)  'tl 

=  '  , 

ft 

WRITE 

(8, 

*) 

rt 

il 

tl,  ' 

.  / 

/ 

WRITE 

(8, 

*) 

't_fin 

=  ',  t_ 

fin,  ';' 

WRITE 

(8, 

*) 

'  total. 

.cpu  =  ' 

,  cputime,  ' 

WRITE 

(8, 

33 

)  'total_cpu  = 

'  ,  cputime. 

11  FORMAT {A4,  E19.7,  Al) 

22  FORMAT (A7,  F16.8,  Al) 

33  FORMAT (All,  F12 . 3 ,  Al) 


!  !  Maximum  error : 

I  I 

value  =  ODO 
error  =  ODO 
DO  ix  =  0,  nx 
xi  =  XX ( ix ) 

DO  jy  =  0,  ny 

yj  =  yy(jy) 

DO  kz  =  0,  nz 
zk  =  zz (kz) 

V  =  func_u { t_f in,  xi,yj,zk) 
value  =  DMAXK  value,  DABS(v)  ) 

error  =  DMAXl (  error,  DABS (  v  -  u(ix,jy,kz)  )  ) 

END  DO 
END  DO 
END  DO 

PRINT  * ,  '  max  value  =  ' ,  value 
PRINT  * ,  '  max  error  =  ' ,  error 
WRITE  (8,  *)  'max_value  =  value, 

WRITE  (8,  *)  'max_error  =  ',  error,  ' ; ' 

!  44  FORMAT (All,  F16.9,  Al ) 

!  55  FORMAT (All,  E16.6,  Al ) 

CLOSE  (UNIT=8) 


CONTAINS 
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!  The  initial  condition: 

I 

FUNCTION  func^uO (x,y, z) 

REAL* 8  func_uO,  x,  y,  z,  temp 
temp  =  DSIN(PI  *  x) 
temp  =  temp  *  DSIN(PI  *  y) 
func_uO  =  temp  *  DSIN(PI  *  z) 

END  FUNCTION  func^uO 

END  PROGRAM  adi3d 

SUBROUTINE  solution ( t_f in,  u) 

!!!!!!!!!!!!  1  !!!!!!  i  !!!!!!!!!!!  i  !!!!!!  I  !!!!!!!!!!!!  i  !  I  !!!!  ! 

New  ADI  solver  for  3D  convection-diffusion  equations 
May  9 ,  1997 

!!!!!!!!!!!!!!!  1  1  !!!!  1  !!!!!!!!!!!!!!  I  !!!  I  !!!!!!!!!!!!!!!!!  ! 

USE  init_data2 
USE  coeffunc 
INTEGER  ktime,  Kmax 
REAL *8  t_fin,  time 
REAL*8  u{0:nx,  0:ny,  0:nz) 

REAL*8  ux(nxl) ,  uy(nyl) ,  uz{nzl) 

REAL* 8  coefOa (dimen,  nxl,  nyl,nzl) 

REAL*8  coef Ob (dimen,  nxl,  nyl,nzl) 

REAL*8  coef Oc (dimen,  nxl,  nyl,nzl) 

REAL* 8  coef la (dimen,  nxl,  nyl,nzl) 

REAL* 8  coef lb (dimen,  nxl,  nyl,nzl) 

REAL*8  coef Ic (dimen,  nxl,  nyl,nzl) 

REAL*8  aax(nx2),  bbx(nxl),  ccx(nx2),  ddx(nxl) 

REAL*8  aay(ny2),  bby{nyl),  ccy(ny2),  ddy(nyl) 

REAL*8  aaz(nz2),  bbz(nzl),  ccz(nz2),  ddz(nzl) 

REAL  cpu,  t array (2),  DTIME 
EXTERNAL  DTIME 
EXTERNAL  coef 
INTRINSIC  IDINT 
INTENT (INOUT)  ::  t_fin 
INTENT (OUT)  ::  u 

cpu  =  DTIME (tarray) 

CALL  coef (  coef Oa, coef Ob, coef Oc ,  coef la, coef lb, coef Ic  ) 

DO  ix  =  0,  nx 
xi  =  XX (ix) 

DO  jy  =  0,  ny 

yj  =  yy(jy) 

DO  kz  =  0 ,  nz 
zk  =  zz{kz) 

u(ix,jy,kz)  =  func_u{t0,  xi,yj,zk) 

END  DO 
END  DO 
END  DO 

cpu  =  DTIME (tarray) 

PRINT  *,  '  cpu  for  the  coefficients  and  initialization  =  cpu 

Kmax  =  IDINT (  ( t_f in-tO ) /dt  ) 

t__fin  =  to  +  dt*Kmax 

IF{Kmax  ,NE.  mt)  PRINT  *,  '  Check:  t_final  might  be  changed.' 

DO  ktime  =  1,  Kmax 

time  =  to  +  dt*ktime 

1,  The  explicit  part,  in  alternating  directions: 


DO  ix 


l,nxl 
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& 

Sc. 


xi  =  XX (ix) 

DO  jy  =  l,nyl 

yj  =  yY(jy) 

DO  kz  =  l,nzl 

uzCkz)  =  coef Oa (3 , ix, jy, kz)  * 
+  coef Ob ( 3 , ix, jy, kz )  * 
+  coef Oc ( 3 , ix, jy, kz )  * 

END  DO 

u(ix,  jy,  l:nzl)  =  uz 
u{ix,  jy,  0)  =  func_u(time, 

u(ix,  jy,  nz)  =  func_u  ( time , 
END  DO 
END  DO 


u{ix, jy,kz-l)  & 
u(ix,jy,kz)  Sc 
u  { ix,  jy,  kz+1) 


xi,  yj ,  zz  (0) ) 
xi,  yj,  zz(nz)) 


DO  kz  =  l,nzl 
zk  =  zz{kz) 

DO  ix  =  1 , nxl 
xi  =  XX (ix) 

DO  jy  =  l,nyl 

uy(jy)  =  coefOa {2 , ix, jy, kz)  *  u ( ix, jy-1 , kz )  & 

Sc  +  coef  Ob  (2 ,  ix,  jy,  kz)  *  u(ix,jy,kz)  & 

Sc  +  coef  Oc  {2  ,  ix,  jy,  kz)  *  u  (ix,  jy+1 ,  kz ) 

END  DO 

u{ix,  l:nyl,  kz)  =  uy 
u(ix,  0,  kz)  =  func_u ( time,  xi,  yy(0),  zk) 
u{ix,  ny,kz)  =  func_u(time,  xi,  yy(ny),zk) 

END  DO 
END  DO 

DO  jy  =  l,nyl 

yj  =  yy(jy) 

DO  kz  =  l,nzl 
zk  =  zz(kz) 

DO  ix  =  1 ,  nxl 

ux(ix)  =  coef0a(l, ix, jy,kz)  *  u ( ix-1 , jy, kz )  & 

Sc  +  coef  Ob  ( 1 ,  ix,  jy,  kz )  *  u(ix,jy,kz)  & 

Sc  +  coefOc  (1,  ix,  jy,  kz)  *  u  ( ix+1 ,  jy,  kz) 

END  DO 

u(l:nxl,  jy,  kz)  =  ux 
u{0,  jy,  kz)  =  func_u{time,  xx{0),  yj ,  zk) 
u(nx,jy,  kz)  =  func_u(time,  xx(nx),yj,  zk) 

END  DO 
END  DO 

u  =  u  +  f(ktime,  :,:,:) 

+  f(ktime,  ix,jy,kz) 


2.  The  implicit  part,  in  alternating  directions: 


DO  j  y  =  1 ,  ny  1 
DO  kz  =  l,nzl 

aax  =  coeflad,  2:nxl,jy,kz) 
bbx  =  coeflbd,  l:nxl,jy,kz) 
ccx  -  coeflcd,  l:nx2,jy,kz) 
ddx  =  u{l:nxl,  jy,  kz) 

ddxd)  =  ddxd)  -  coeflad,  1,  jy,  kz)  *u  (0  ,  jy,  kz) 
ddx(nxl)  =  ddx(nxl)  -  coeflcd,  nxl ,  jy,  kz)  *u  (nx,  jy,  kz) 
CALL  tridiag(nxl,  aax,  bbx,  ccx,  ddx,  ux) 

IF(ktime==l  .AND.  jy==l  .AND.  kz==l)  & 

Sc  PRINT  *,  '  after  calling  tridiag 

ud:nxl,  jy,  kz)  =  ux 
END  DO 
END  DO 
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!  IF{ktime  .EQ.  1)  PRINT  '  one  implicit  step  survived 

DO  kz  =  l,nzl 
DO  ix  =  l,nxl 

aay  =  coefla(2,  ix,2:nyl,kz) 
bby  =  coeflb{2,  ix,l:nyl,kz) 
ccy  =  coeflc(2,  ix,l:ny2,kz) 
ddy  =  u{ix,  l:nyl,  kz) 

ddy(l)  =  ddy(l)  -  coefla{2,  ix, 1 , kz ) *u ( ix, 0 , kz ) 
ddy(nyl)  =  ddy{nyl)  ~  coeflc(2,  ix, nyl , kz) *u ( ix, ny, kz ) 
CALL  tridiag(nyl,  aay,  bby,  ccy,  ddy,  uy) 
u(ix,  l:nyl,  kz)  =  uy 
END  DO 
END  DO 

DO  ix  =  1 , nxl 
DO  jy  =  l,nyl 

aaz  =  coeflaO,  ix,jy,2:nzl) 
bbz  =  coeflbO,  ix,jy,l:nzl) 
ccz  =  coeflcO,  ix,jy,l:nz2) 
ddz  =  u(ix,  jy,  Irnyl) 

ddz(l)  =  ddz{l)  -  coeflaO,  ix,  jy ,  1 )  *u  (ix,  jy,  0 ) 
ddz(nzl)  =  ddz(nzl)  -  coeflcO,  ix,  jy,  nzl)  *u  (ix,  jy,nz) 
CALL  tridiag(nzl,  aaz,  bbz,  ccz,  ddz,  uz) 
u{ix,  jy,  l:nzl)  =  uz 
END  DO 
END  DO 

WHERE  (u  <  O.ODO)  u  =  0 . ODO 

WHERE  (u{:,:,:)  <  O.ODO)  u(:,  :,  :)  =  O.ODO 

END  DO 
RETURN 

END  SUBROUTINE  solution 

SUBROUTINE  coef (  coef Oa , coef Ob , coef Oc ,  coef la, coef lb, coef Ic  ) 
SUBROUTINE  coef (  coef Oa, coef Ob, coef Oc,  coef la, coef lb, coef Ic ,  f  ) 
!!!!  1  !!!!!!!!  i  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I  !!!!!!!  ! 

!  The  coefficients  for  solving  the  ID  sub-problems 

!  May  9 ,  1997 

I  j  I  I  I  I  I  I  I  1  !  j  I  I  1  !  !  I  I  i  I  I  I  I  I  I  I  I  I  !  I  I  I  j  I  I  I  !  I  I  !  !  !  I  !  !  I  !  !  !  I  I  !  I  !  I  I  I  !  I 

USE  init_data2 
USE  coeffunc 

REAL*8  coef Oa (dimen,  nxl,  nyl,nzl) 

REAL *8  coef Ob (dimen,  nxl,  nyl,nzl) 

REAL*8  coef Oc (dimen,  nxl,  nyl,nzl) 

REAL* 8  coef la (dimen,  nxl,  nyl,nzl) 

REAL *8  coef lb (dimen,  nxl,  nyl,nzl) 

REAL *8  coef Ic (dimen,  nxl,  nyl,nzl) 

REAL* 8  f(0:nx,  0:ny,  0:nz) 

REAL*8  dxyz (dimen),  dneg ( dimen ) ,  dpos (dimen) 

REAL*8  dxyz_neg,  temp,  temp_c,  vect_b (dimen) 

LOGICAL  judge (dimen) 

COMMON  /coef_right/  convection,  diffusion,  dxyz 
INTRINSIC  DABS,  DSIN,  DCOS 

INTENT(OUT)  ::  coef Oa , coef Ob, coef Oc ,  coef la, coef lb, coef Ic 
DO  i  =  1 ,  dimen 

dxyz(i)  =  dif fus (i) *dtx2 (i) /dx(i) 
temp  =  2  *  dxyz(i) 
dneg(i)  =  1  -  temp 
dpos(i)  =  1  +  temp 


END  DO 
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DO  ix 


DO  jy 

yj 

DO  kz 


1 ,  nxl 
XX ( ix) 
1  /  nyl 

yy  ( jy) 

1 ,  nzl 
zz (kz ) 


vect_b  =  func_b{  xi ,  yj ,  zk  ) 

temp_c  =  func_c(  xi,  yj,  zk  )  *  dtdim2 

judge  =  (vect_b  .LE.  ODO) 


DO  i  =  1 ,  dimen 

temp  =  vect_b(i)  *  dtx2(i) 
temp_b  =  DABS  (temp) 
dxyz_neg  =  -  dxyz ( i ) 

IF  ( judge (i))  THEN 

coefOa(i,  ix,jy,kz)  =  dxyz(i) 

coefOc(i,  ix,jy,kz)  =  dxyz(i)  +  temp 

coefla(i,  ix,jy,kz)  =  dxyz_neg  +  temp 
coeflc(i,  ix,jy,kz)  =  dxyz_neg 
ELSE 

coefOa(i,  ix,jy,kz)  =  dxyz(i)  -  temp 

coefOc(i,  ix;jy,kz)  =  dxyz(i) 

coefla(i,  ix,jy,kz)  =  dxyz_neg 

coeflc(i,  ix,jy,kz)  =  dxyz_neg  -  temp 
END  IF 

coefOb(i,  ix,jy,kz)  =  dneg(i)  +  temp_b  +  temp_c 

coeflb(i,  ix,jy,kz)  =  dpos(i)  +  temp_b  -  temp_c 

END  do' 


END  DO 
END  DO 
END  DO 

RETURN 

END  SUBROUTINE  coef 

SUBROUT INE  f _add ( w ,  j  j  udge ,  dd_ j  udge , dd ,  f  f ) 
REAL  *8  w,  dd__judge,  dd,  ff 
LOGICAL  j  j  udge 

INTENT  (IN)  ::  w,  j  judge,  dd_  judge,  dd 
INTENT (INOUT)  ::  ff 
IF(w  .NE.  ODO)  THEN 
IF (j judge)  THEN 

ff  =  ff  +  w  *  dd_ judge 


ELSE 

ff  =  ff  +  w  * 
END  IF 
END  IF 
RETURN 

END  SUBROUTINE  f  add 


SUBROUTINE  tridiag (N, A, B, C , D, X) 
INTEGER  N 

REAL*8  A(N-l) ,B(N) ,C(N-1) ,D(N) ,X(N) 
INTENT (IN)  ::  N,  A,C 
INTENT (INOUT)  ::  B,D,  X 
DO  2  I  =  2,N 

XMULT  =  A(I-l) /B(I-l) 

B(I)  =  B(I)  -  XMULT*C(I-1) 

D{I)  =  D(I)  -  XMULT*D(I-1) 

CONTINUE 
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X{N)  =  D{N)/B(N) 

DO  3  I  =  N-1,1,-1 

X(I)  =  (D{I)  -  C(I) *X(I+1) ) /B(I) 

3  CONTINUE 
RETURN 

END  SUBROUTINE  tridiag 
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!!!!!!!!!!!!!!!!!  I  !!  1  !!!  I  !  i  I  I  1  I  M  I  I  1  I  I  t  I  1  j  I  j  I  I  I  j  I  I  I  I  I  ! 
On-line  computation  of  a  3-D  filtering  problem 
Chuanxia  Rao 
January  1997 
April  1996 

!  I  !  !  !  I  I  I  !  »  1  1  I  !  !  I  1  I  J  I  M  !  !  i  I  I  I  !  !  !  I  I  j  t  I  I  I  I  J  I  I  I  I  I  I  I  I  I  J  I  I  I  t 
Dependencies : 

data_h  --  computed  here,  simple 

data_v  --  off line. data 

data_z  --  observed. data 

1  I  1  !  !  I  I  !  !  I  !  !  I  I  I  !  !  !  !  I  !  !  I  I  I  1  1  I  I  I  I  [  [  1  1  1  I  I  I  I  I  I  M  I  I  I  I  I  [  I  1  I  [ 


Modules  for  initial  data: 

MODULE  init_dim 
INTEGER  dimen ,  dim_ob 
PARAMETER  (dimen  =  3) 

PARAMETER  (dim_ob  =  2) 

END  MODULE  init^dim 

MODULE  init_datal 
USE  init_dim 

REAL*8,  PARAMETER  : :  PI  =  3.14159265358979323 
REAL *8  xi,  yj ,  zk 
INTEGER  ix,  jy,  kz 

INTEGER,  PARAMETER  ::  nx=32,  ny=32,  nz=32 
INTEGER,  PARAMETER  ::  nx=64,  ny=64,  nz=64 
INTEGER,  PARAMETER  : :  nxl=nx-l,  nyl=ny-l,  nzl=nz-l 
INTEGER,  PARAMETER  ::  nx2=nx-2,  ny2=ny-2,  nz2=nz-2 
INTEGER,  PARAMETER  : :  mt  =  100 
INTEGER,  PARAMETER  : :  mt  =  200 
INTEGER,  PARAMETER  : :  msubt  =  4 
END  MODULE  init_datal 

MODULE  init_data2 

USE  init_datal 

REAL  *  8  di f  f us ( dimen ) 

REAL *8  dt,  dx (dimen) 

REAL*8  dtdim2,  dtx2 (dimen) 

REAL* 8  to,  tl 

REAL* 8  xa,  xb,  ya,  yb,  za,  zb 
REAL* 8  xx(0:nx) 

REAL* 8  yy(0:ny) 

REAL*8  zz(0:nz) 

PARAMETER  (t0=0D0,  tl=lD0) 

PARAMETER  (xa=-1.0D0,  xb=:1.0D0) 

PARAMETER  (ya=-1.0D0,  yb=1.0D0) 

PARAMETER  (za=-1.0D0,  zb=1.0D0) 

PARAMETER  (  dx  =  ( /  (xb-xa) /nx,  (yb-ya) /ny,  (zb-za)/nz  /)  ) 

PARAMETER  (dt  =  IDO/mt) 

PARAMETER  (dt  =  (tl-t0)/mt) 

PARAMETER  (dtdim2  =  dt/dimen/2) 

PARAMETER  (  dtx2  =  dt/dx/2  ) 

PARAMETER  (subdt  =  dt /msubt)  ! ! (subdt  =  dt) 

PARAMETER  (  xx  =  xa  +  dx(l)  *  (/(i,  i=0,nx)/)  ) 

PARAMETER  (  yy  =  ya  +  dx(2)  *  (/(j,  j=0,ny)/)  ) 

PARAMETER  (  zz  =  za  +  dx(3)  *  (/(k,  k=0,nz)/)  ) 

PARAMETER  (  diffus  =  (/  1D0/6D0,  1D0/6D0,  1D0/6D0  /)  ) 

PARAMETER  (  diffus  =  (/  IDO,  IDO,  IDO  / )  ) 

DATA  diffus  /IDO,  2D0,  3D0/  !!  diffusion  parameters 
END  MODULE  init_data2 
MODULE  init_noise 
USE  init_dim 


J 
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REALMS  delta (dimen) 

REAL *8  sigma {dim_ob) 

DATA  sigma  /2.3D-1,  3.0D'-2/  !!  noise  parameters  in  observation 

DATA  delta  /1.2D-1,  2.0D~2,  l.OD-2/  !!  noise  parameters  in  signal 
END  MODULE  init_noise 

MODULE  init_dens 
USE  init_dim 

REAL*8  ■  xmean (dimen) ,  var_inv (dimen, dimen) 

! !  initial  mean  and  variance  inverse 
DATA  xmean  /2.3D-1,2,  3.0D~2/ 

DATA  var_inv  /1.2D1,0,0,  0,2.0D2,0,  0,0,1.0D2/ 

END  MODULE  init_dens 

MODULE  ini t_onl ine 

CHARACTER  (*)  ::  file_z,  fmt_z,  status_z 

CHARACTER  (*)  ::  file_dens,  fmt_dens,  status_dens 

INTEGER  Mmax 
PARAMETER  (Mmax  =  50) 

PARAMETER  (file_z  =  'data_z') 

PARAMETER  ( fmt_z  =  ' ( 2F17 . 9 ) ' ) 

PARAMETER  (status_z  =  'OLD') 

PARAMETER  (file_dens  =  'density.m' ) 

PARAMETER  ( f mt_dens  =  ' (2X,  E22.14)') 

PARAMETER  (status_dens  =  'REPLACE') 

END  MODULE  init_online 

MODULE  init_out 

CHARACTER  (*)  file_out  !,  fmt_out,  status_out 
PARAMETER  (file_out  =  ' result . temp ' ) 

END  MODULE  init_out 

MODULE  coeffunc 
USE  init_datal 
CONTAINS 

The  coefficient  functions  for  convection  (or  drift)  term: 

FUNCTION  func_.b(x,y,  z) 

REAL  *  8  f unc_b ( dimen ) 

REAL* 8  X,  y,  z 
func_b(l)  =  -  DSIN(PI  *  x) 
func_b(2)  =  -  DSIN(PI  *  y) 
func_b(3)  =  ~  DSIN(PI  *  z) 

END  FUNCTION  func_b 

The  coefficient  function  for  u  (zero-th  order  term): 

FUNCTION  func_c(x,y,z) 

REAL*8  func_c,  x,  y,  z,  temp 
temp  =  DCOS(PI  *  x) 
temp  =  temp  +  DCOS(PI  *  y) 

temp  =  temp  +  DCOS(PI  *  z) 

temp  =  temp  -  PI 

func_c  =  PI  *  temp 

END  FUNCTION  func_c 

The  force  function  f: 

FUNCTION  func_f (t,x,y, z) 

REAL* 8  func_f,  t,  x,  y,  z,  temp 

END  FUNCTION  func_f 


FUNCTION  func_u(t,  x,y,z) 
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REAL* 8  func_u,  t,  x,  y,  z,  temp 
REAL*8,  PARAMETER  : ;  PI  =  3.14159265358979323 
INTRINSIC  DSIN,  DEXP 
temp  =  DSIN(PI  *  x) 
temp  =  temp  *  DSIN(PI  *  y) 
temp  =  temp  *  DSIN(PI  *  z) 
func_u  =  temp  *  DEXP ( -3D0/2 *PI*PI* t ) 
func_u  -  temp  *  DEXP {-4*PI*PI*t ) 

END  FUNCTION  func_u 
END  MODULE  coeffunc 

FUNCTION  func_h{x) 

REAL  *  8  f  unc_h  ( dim__ob ) 

REAL *8  X (dimen) 

REAL* 8  temp,  tempi,  temp2 
INTRINSIC  ASIN,  ACOS,  SQRT 
temp  =  x{l) **2  +  x(2) **2 
tempi  =  x{l)  /  SQRT (temp) 
temp2  =  x(3)  /  SQRT ( temp+x (3 ) **2 ) 
f unc_h ( 1 )  =  ACOS ( tempi ) 
func_h(2)  =  ASIN{temp2) 

END  FUNCTION  func_h 

PROGRAM  filter3adi 

Main  program  begins : 

USE  init_cpu 

USE  init_data2 

USE  init_dens 

USE  init_noise 

USE  init_online 

REAL *8  Z(dim_ob,  Mmax) 

REAL *8  hh(dim_ob,  0 : nx, 0 : ny , 0 : nz ) 

REAL*8  p(0:l,  0 : nx, 0 : ny , 0 : nz ) 

REAL* 8  u ( 0 :nx, 0 : ny, 0 :nz ,  -Lx : Lx, -Ly : Ly, ~Lz : Lz) 

REAL  *  8  pmax ,  pmin 
REAL *8  tempi,  temp,  tau 

REAL*8  tempi jk (dim_ob) ,  ob_var_inv (dim_ob, dim_ob) 
INTRINSIC  ABS,  EXP,  MAX,  MIN  !!  DABS,  DEXP,  DMAXl 
INTRINSIC  MATMUL,  DOT_PRODUCT 

PRINT  *,  'Reading  data  and  initialization' 

PRINT  *,  'and  off-line  calculations,  please  wait  ...' 

OPEN  {UNIT=8,  FILE=file_z,  STATUS=status_z ) 

READ  (UNIT=8,  FMT=fmt_z)  {  ( Z ( i , m) , i=l , dim_ob) ,  m=l,Mmax) 

CLOSE  (UNIT=8) 

p(0,  :,:,:)  =  func_pO ( { /xx , yy , zz/ ) ) 

DO  i  =  0 ,  nx 
DO  j  =0,  ny 
DO  k  =  0,  nz 

hh{:,  i,j,k)  =  func_h(  (/xx(i),  yy{j),  zz(k)/)  ) 
p(0,  i,j,k)  =  func_p0(  (/xx(i),  yy(j),  zz(k)/)  ) 
END  DO 
END  DO 
END  DO 

PRINT  *,  'On-line  calculations  now  . . . ' 

PRINT  * ,  ' Mmax  = ' ,  Mmax 

tau  =  40 

tau  =  14.14 
tau  =  21.21 


exp(-40)  =  4.25e-18 
!  =  sqrt(2*100) 

!  =  sqrt(2*225) 
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ob_var_inv  =  0 
DO  ii  =  1,  diin_ob 

temp  =  IDO  /  sigma (ii) 
ob_var_inv(ii, ii)  =  temp  *  temp 
END  DO 

time_loop:  DO  m  =  0,  Mmax-1 
ml  =  m  +  1 

CALL  of f_line {xx,yy, zz ,  indx, indy , indz ,  u) 
pmax  =  0.0 
pmin  =  0.0 
p(l,  =  0.0 

p  (ml,  =  0.0 

DO  ix  =  0,nx 
DO  jy  =  0,ny 
DO  kz  =  0,nz 

tempijk  =  Z(:,ml)  -  hh ( : , jx, jy, j z) 

tempi  =  DOT_PRODUCT ( tempi jk,  MATMUL (ob_var_inv, tempi jk) /2 ) 

IF (tempi  .LE.  tau)  THEN 

tempi  =  EXP (“tempi) 

temp  =  u(ix,jy,kz) 

p(ml,  jx,jy,jz)  =  tempi  *  temp 

p(l,  jx,jy,jz)  =  tempi  *  temp 

pmax  =  MAX  (pmax,  p(l,  jx,jy,jz)) 

pmin  =  MIN  (pmin,  p(l,  jx,jy,jz)) 

END  IF 

END  DO 
END  DO 
END  DO 

IF (pmin  .LT.  0.0)  THEN 

PRINT  * ,  '  p_min  <0  ...  ' 

p(l,  ;,:,:)  =p(l,  :,:,:)  -  pmin 
pmax  =  pmax  -  pmin 
END  IF 

IF (pmax  .GT.  0.0)  THEN 

WHERE (p(l, :,:,:)>0.0)  & 

&  p(l,  :,:,:)  =  p(l,  :,:,:)  /  pmax 

ELSEWHERE 
END  WHERE 
END  IF 

IF (ml  .LT.  Mmax)  p(0,  :,:,:)  =p(l,  :,:,:) 

END  DO  time_loop 

PRINT  *,  'Saving  results  of  the  last  step 

OPEN  (UNIT=9,  FILE=file_dens,  STATUS=s tatus_dens ) 

WRITE  (UNIT=9,  FMT='(5A)')  'p  =  [ ' 

DO  i  =  0,  nx 
DO  j  =0,  ny 
DO  k  =  0,  nz 

temp  =  p (1,  i, j , k) 

IF (temp  .GE.  0.99D-99)  THEN 

WRITE  (UNIT=9,  FMT=fmt_dens )  temp 
ELSE 

WRITE  (UNIT=9,  FMT= ' ( 4X, A3 ) ' )  '0.0' 

END  IF 
END  DO 
END  DO 
END  DO 

&  (((p(Mmax,  i,j,k),  i=0,nx),  j=0,ny),  k=0,nz) 

WRITE  (UNIT=9,  FMT='(2A)')  ' ] ; ' 

CLOSE  (UNIT=9) 


k=0 , nz) 


CONTAINS 


The  initial  density  function  pO  is  defined  here: 
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FUNCTION  func_pO{x) 

REAL  *8  func^pO,  t  emp 

REAL*8,  DIMENSION  (dimen)  ::  x,  tempO ,  tempi 
INTRINSIC  MATMUL,  DOT_PRODUCT,  EXP 
tempo  =  X  -  xmean 
tempi  =  MATMUL (var_inv,  tempO ) 
temp  =  MATMUL {TRANSPOSE ( tempo ), tempi) 
temp  =  DOT_PRODUCT ( tempo ,  tempi) 
func_pO  =  EXP (-temp) 

END  FUNCTION  func_j)0 

END  PROGRAM  filterSadi 


SUBROUTINE  adi3d 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I  !!!!!!!!!!!  1  !  1  !!!!!!!  ! 

New  ADI  for  convection-diffusion  equations 
Chuanxia  Rao 
May  1997 
Feb  1996 

!!!!!!!!!!!!!!!!!!  I  !!!!!!!  I  !!!!  I  !!!!!!!!!!!!!!!  I  !!!!!  ! 
SUBROUTINE  solution ( t__f in,  u) 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  I  !!!!!!!!!!  1  1  !!!  1 
New  ADI  solver  for  3D  convection-diffusion  equations 
May  9,  1997 

I  I  I  I  M  I  I  I  I  I  I  ]  I  t  I  I  I  I  I  I  I  I  I  I  [  I  I  I  I  I  ]  I  I  I  I  I  I  I  I  I  I  I  I  I  I  M  M  I  I  I  I  I  I  I  I  I 


USE  init_data2 
USE  coeffunc 
INTEGER  ktime,  Kmax 
REAL *8  t_fin,  time 
REALMS  u(0:nx,  0:ny,  0;nz) 

REAL*8  ux(nxl),  uy(nyl),  uz(nzl) 
REAL*8  coefOa (dimen,  nxl,  nyl,nzl) 
REAL* 8  coef Ob (dimen,  nxl,  nyl,nzl) 
REAL* 8  coef Oc (dimen,  nxl,  nyl,nzl) 
REAL* 8  coef la (dimen,  nxl,  nyl,nzl) 
REAL* 8  coef lb (dimen,  nxl,  nyl,nzl) 
REAL* 8  coef Ic (dimen,  nxl,  nyl,nzl) 
REAL*8  aax(nx2),  bbx(nxl),  ccx(nx2), 
REAL* 8  aay(ny2),  bby(nyl) ,  ccy(ny2), 
REAL*8  aaz(nz2),  bbz{nzl),  ccz(nz2), 
REAL  cpu,  tarray(2),  DTIME 
EXTERNAL  DTIME 
EXTERNAL  coef 
INTRINSIC  IDINT 
INTENT (INOUT)  ::  t_fin,  u 
INTENT (OUT)  u 


ddx(nxl) 
ddy  (nyl) 
ddz (nzl) 


cpu  =  DTIME (tarray) 

CALL  coef (  coef Oa , coef Ob, coef Oc ,  coef la , coef lb, coef Ic  ) 

DO  ix  =  0,  nx 
xi  =  XX (ix) 

DO  jy  =  0,  ny 

yj  =  yy{Dy) 

DO  kz  =  0 ,  nz 
zk  =  zz{kz) 

u{ix,jy,kz)  =  func_u(t0,  xi,yj,zk) 

END  DO 
END  DO 
END  DO 

cpu  =  DTIME (tarray) 

PRINT  *,  '  cpu  for  the  coefficients  and  initialization  =  cpu 
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Kmax  =  IDINT(  { t_f in-tO ) /dt  ) 
t_fin  =  to  +  dt*Kmax 

IF (Kmax  .NE.  rat)  PRINT  '  Check:  t_final  might  be  changed.' 

DO  ktirae  =  1,  Kmax 

time  =  to  +  dt*ktime 

!  t 

!!  1.  The  explicit  part,  in  alternating  directions: 

{  I 

DO  ix  =  1 , nxl 
xi  =  XX (ix) 

DO  jy  =  l,nyl 

yj  =  yy(jy) 

DO  kz  =  l,nzl 

uz{kz)  =  coef0a(3,ix, jy,kz)  *  u (ix, jy, kz-1)  & 

&  +  coef0b{3 , ix, jy, kz)  *  u(ix,jy,kz)  & 

Sc  +  coefOc  (3,  ix,  jy,kz)  *  u  { ix,  jy,  kz  +  1 ) 

END  DO 

u(ix,  jy,  l:nzl)  =  uz 

u{ix,  jy,  0)  =  func_u(time,  xi,  yj ,  zz(0)) 

u(ix,  jy,  nz)  =  func_u{time,  xi,  y j ,  zz{nz)) 

END  DO 
END  DO 

DO  kz  =  l,nzl 
zk  =  zz{kz) 

DO  ix  =  1 , nxl 
xi  =  XX (ix) 

DO  jy  =  l,nyl 

uy(jy)  =  coef0a(2, ix, jy,kz)  *  u (ix, jy-l, kz)  & 

+  coef Ob (2 , ix, jy, kz)  *  u(ix,jy,kz)  & 

+  coef  Oc  (2  ,  ix,  jy,  kz)  *  u  ( ix,  jy-fl ,  kz) 

END  DO 

u ( ix ,  1 : nyl ,  kz )  =  uy 
u(ix,  0,  kz)  =  func_u(time,  xi,  yy(0) ,  zk) 
u{ix,  ny,kz)  =  func_u(time,  xi,  yy(ny),zk) 

END  DO 
END  DO 

DO  jy  =  1/nyl 

yj  =  yy(jy) 

DO  kz  =  l,nzl 
zk  =  zz (kz ) 

DO  ix  =  1 , nxl 

ux(ix)  =  coefOad,  ix,  jy,kz)  *  u  ( ix-1 ,  jy,  kz )  & 

+  coef Ob (1 , ix, jy, kz )  *  u(ix,jy,kz)  & 

+  coef Oc ( 1 , ix, jy, kz )  *  u(ix+l, jy,kz) 

END  DO 

u(l:nxl,  jy,  kz)  =  ux 
u(0,  jy,  kz)  =  func_u(time,  xx(0) ,  yj ,  zk) 
u(nx,jy,  kz)  =  func_u(time,  xx(nx),yj,  zk) 

END  DO 
END  DO 

u  =  u  +  f(ktime,  :,:,:) 

+  f(ktime,  ix,jy,kz) 


!!  2.  The  implicit  part,  in  alternating  directions: 


DO  jy  =  l,nyl 
DO  kz  =  l,nzl 

aax  =  coeflad,  2:nxl,jy,kz) 
bbx  =  coeflb(l,  l:nxl,jy,kz) 
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ccx  =  coeflc(l,  l:nx2,jy,kz) 
ddx  =  u(l:nxl,  jy,  kz) 

ddx(l)  =  ddx(l)  -  coeflad,  1 ,  jy ,  kz )  ( 0  ,  jy,  kz ) 

ddx{nxl)  =  ddx{nxl)  -  coeflcd,  nxl ,  jy ,  kz )  *u  (nxoV/ kz ) 
CALL  tridiag(nxl,  aax,  bbx,  ccx,  ddx,  ux) 
ud:nxl,  jy,  kz)  =  ux 
END  DO 
END  DO 

IF(ktime  .EQ.  1)  PRINT  *,  '  one  implicit  step  survived 

DO  kz  =  l,nzl 
DO  ix  =  1 , nxl 

aay  =  coefla(2,  ix,2:nyl,kz) 
bby  =  coeflb(2,  ix,l:nyl,kz) 
ccy  =  coeflc(2,  ix,l:ny2,kz) 
ddy  =  u { ix ,  1 : nyl ,  kz ) 

ddy(l)  =  ddy(l)  ~  coefla{2,  ix, 1 , kz ) *u ( ix, 0 , kz ) 
ddy (nyl)  =  ddy (nyl)  -  coeflc(2,  ix,nyl, kz) *u {ix,ny, kz) 
CALL  tridiag{nyl,  aay,  bby,  ccy,  ddy,  uy) 
u ( ix ,  1 : nyl ,  kz )  =  uy 
END  DO 
END  DO 

DO  ix  =  l,nxl 
DO  jy  =  l,nyl 

aaz  =  coeflaO,  ix,jy,2:nzl) 
bbz  =  coeflb(3,  ix,jy,l:nzl) 
ccz  =  coeflcO,  ix,jy,l:nz2) 
ddz  =  u(ix,  jy,  l:nyl) 

ddzd)  =  ddzd)  -  coefla{3,  ix,  jy,  1)  *u  (ix,  jy,  0) 
ddz(nzl)  =  ddz{nzl)  -  coeflc(3,  ix, jy,nzl) *u(ix, jy,nz) 
CALL  tridiag(nzl,  aaz,  bbz,  ccz,  ddz,  uz) 
u(ix,  jy,  l:nzl)  =  uz 
END  DO 
END  DO 

WHERE  (u  <  O.ODO)  u  =  0 . ODO 

WHERE  (u(:,:,:)  <  O.ODO)  u(:,  :,  :)  =  0 . ODO 

END  DO 
RETURN 

END  SUBROUTINE  adi3d 
END  SUBROUTINE  solution 


SUBROUTINE  coef (  coef Oa, coef Ob, coef Oc ,  coef la, coef lb, coef Ic  ) 
SUBROUTINE  coef (  coef Oa, coef Ob, coef Oc ,  coef la, coef lb, coef Ic , 
1  !!!  I  !!!!!!!!  1  !!!!!!!!!!!!!!!!!!!!  1  !!!!!!  I  !!!!!!!!!!!!  I  !!!  ! 

The  coefficients  for  solving  the  ID  sub-problems 

May  9,  1997 
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I 


USE  init_data2 
USE  coeffunc 
REAL*8  coef Oa (dimen, 
coef Ob (dimen. 


REAL *8 
REAL *8 
REALMS 
REAL *8 
REAL *8 
REAL* 8 
REAL* 8 
REAL *8 


nxl,  nyl,nzl) 
nxl,  nyl,nzl) 


coef Oc (dimen,  nxl,  nyl,nzl) 
coef la (dimen,  nxl,  nyl,nzl) 


nxl,  nyl,nzl) 
nxl,  nyl,nzl) 


coef lb (dimen, 
coef Ic (dimen, 
f(0:nx,  0;ny,  0:nz) 
dxyz (dimen),  dneg( dimen),  dpos (dimen) 
dxyz_neg,  temp,  temp_c,  vect_b (dimen) 


LOGICAL  j  udge ( dimen ) 


f  ) 
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COMMON  /coef_right/  convection,  diffusion,  dxyz 
INTRINSIC  DABS,  DSIN,  DCOS 

INTENT(OUT)  ::  coef Oa , coef Ob, coef Oc ,  coef la, coef lb, coef Ic 


DO  i  =  1,  dimen 

dxyz(i)  =  dif fus (i) *dtx2 (i) /dx(i) 
temp  =  2  *  dxyz(i) 
dneg ( i )  =  1  -  temp 
dpos(i)  =  1  +  temp 
END  DO 


DO  ix  =  1,  nxl 
xi  =  XX (ix) 
DO  jy  =  1,  nyl 

yj  =  yy(jy) 

DO  kz  =  1,  nzl 
zk  =  zz (kz ) 


vect_b  =  func_b(  xi ,  yj ,  zk  ) 

temp_c  =  f unc_c  (  xi ,  yj  ,  zk  )  dtdim2 

judge  =  (vect_b  .LE.  ODO) 


DO  i  =  1 ,  dimen 

temp  =  vect_b(i)  *  dtx2(i) 
temp_b  =  DABS (temp) 
dxyz__neg  =  -  dxyz  ( i ) 

IF  (judge(i))  THEN 

coefOa(i,  ix,jy,kz)  =  dxyz(i) 

coefOc(i,  ix,jy,kz)  =  dxyz(i)  +  temp 

coefla(i,  ix,jy,kz)  =  dxyz_neg  +  temp 

coeflc(i,  ix,jy,kz)  =  dxyz_neg 

ELSE 

coefOa(i,  ix,jy,kz)  =  dxyz(i)  -  temp 

coefOc(i,  ix,jy,kz)  =  dxyz(i) 

coefla(i,  ix,jy,kz)  =  dxyz_neg 

coeflc(i,  ix,jy,kz)  =  dxyz_neg  -  temp 

END  IF 

coefOb(i,  ix,jy,kz)  =  dneg(i)  +  temp_b  +  temp_c 

coeflb(i,  ix,jy,kz)  =  dpos(i)  +  temp_b  -  temp_c 

END  DO 


END  DO 
END  DO 
END  DO 


RETURN 

END  SUBROUTINE  coef 

SUBROUTINE  f_add(w,  j judge,  dd_judge,dd,  ff) 
REAL* 8  w,  dd_judge,  dd,  ff 
LOGICAL  j  j  udge 

INTENT  (IN)  w,  j  judge,  dd_ judge,  dd 
INTENT (INOUT)  ::  ff 
IF(w  .NE.  ODO)  THEN 
IF (j judge)  THEN 

ff=ff+w*  dd_ judge 
ELSE 

ff  =  ff  +  w  *  dd 
END  IF 
END  IF 
RETURN 

END  SUBROUTINE  f_add 
SUBROUTINE  tridiag (N, A, B , C , D, X) 
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INTEGER  N 

REAL*8  A{N-1) , B (N) ,C(N-1) ,D(N) ,X{N) 
INTENT (IN)  N,  A,C 
INTENT (INOUT)  ::  B,D,  X 
DO  2  I  =  2,N 

XMULT  =  A(I-l) /B(I-l) 

B{I)  =  B(I)  -  XMULT*C(I-1) 

D(I)  =  D(I)  -  XMULT*D(I-1) 

2  CONTINUE 

X(N)  =  D(N) /B(N) 

DO  3  I  =  N-1,1,-1 

X(I)  =  (D(I)  -  C(I) *X(I+1) ) /B(I) 

3  CONTINUE 
RETURN 

END  SUBROUTINE  tridiag 


